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ABSTRACT 

We present the first numerical result of fully general relativistic axisymmetric simulations for the 
collapse of a rotating high-entropy stellar core to a black hole and an accretion disk. The simulations 
are performed taking into account the relevant microphysics. We adopt as initial condition a spherical 
core with constant electron fraction (Y e = 0.5) and entropy per baryon s — 8 fcg, and angular 
velocity is superimposed. In the early phase, the core collapses in a homologous manner. Then, it 
experiences a weak bounce due to the gas pressure of free nucleons. Because the bounce is weak, the 
core collapses eventually to a black hole. Subsequent evolution depends on initial angular velocity. 
When the rotation is not fast, a geometrically thin (but optically thick) accretion disk is formed, 
and shock waves are formed in the inner part of the disk. For the moderately rotating case, the thin 
accretion disk expands eventually to be a geometrically thick torus after sufficient accumulation of the 
thermal energy generated at the shocks. Furthermore, convection occurs inside the torus. Neutrino 
luminosities vary violently with time because of the convective motion. For the rapidly rotating case, 
by contrast, a geometrically thick torus is formed soon after the black hole formation, and convective 
activity is weak due to the presence of epicyclic mode. 

Subject headings: black hole physics - gamma rays:bursts - accretion, accretion disks - stars: rotation 



1. INTRODUCTION 

Gamma-ray bursts (GRBs) have been one of the most 
outstanding ph enomena in the univers e since their dis- 
covery in 1967 (jKlebesadel et all 1 19731 ) because of their 
huge energy emitted in a short timescale (isotopic equiv- 
alent luminosities of 10 49 -10 52 ergs/s in short duration 
of ~ 0.01-1000 s) and in addition, violent time variabil- 
ity of St ~ 1 ms in time profiles of gamma-ray emis- 
sion. GRBs are basically divided, in terms of their du- 
ration, into short bursts (SGRBs), for which duration is 
shorter than 2 s, and long bursts (LGRBs), for which 
duration is longer than 2 s. Recent observations have 
found GRBs with o verlapped features of t he two popula- 
tions ([Gehrels et al.ll2006HGal-Yam et al.ll2006[ ). and it is 
also suggested that a new classifica tion may be necessary 
(|Zhang et all 120091: iLu et all 12010ft . However, the large 
amount of energy release, short duration, and variability 
timescale indicate that GRBs may be universally associ- 
ated with accreti on processe s onto a compact object of 
stellar- mass size (|PiranllT999l ). Because a rotating black 
hole is the most efficient converter of gravitational bind- 
ing energy in nature, it is now widely believed that many 
of central engines of GRBs are composed of a rotating 
black hole surrounded by a massive and hot accretion 
disk. 

Although progenitors of GRBs have not been fully 
clarified yet, there are accumulating observational ev- 
idences that LGRBs are associated w ith collapse of 
massive stars (Woosle v fe Blooml |2006[). (For re views 
on progenitors of SGRBs, see e.g.. iNakan (|2007l ) and 
ILee fc Ramirez-Ru iz (2007)). The first solid evidence 
for the connection between LGRBs and supernovae 
came from spectroscopic identification of a supernova 
component (SN2003dh) in the afterglow of GRB030329 



dHiorth et al.l l2003t iStanek et ail 120031: iKawabata et al.1 
1200311 . To date, at least six other connections be- 
tween LGRBs and supernovae have been re p orted : 
GRB980425 with S N1998bw (IGalama et al.l H9981 
Kul karni et all [1998); XRF020903 (ISoderber g et al.l 
20051) : GRB021211 with SN200 21t dDella Valle et"a!1 
pool GRB03120 3 with SN20031w (iMalesani et al.ll2004t 
Cobb et al.1 12004 iThomsen etafl 120041: iGal-Yam et al 
I2004D: GRB050525a with SN2005nc dDella Valle et"a! 



2006 rj) and GRB060218 with SN2006ai (iCampana et al 
2006 1 : i Pian et al.ll2006t iMirabal et aTll2006l:lModiaz et al 
l2006t lSollerman et all l2006[): All the GRB-associated su- 
pernovae are Typelb/c. In addition there are a wide 
varie ty of circumstance evidences ([Wooslev fc Blooml 
2006): E.g., observed association of afterglows of 
LGRBs with star forming re g ions in their host 
galaxies (IChristensen et al.l 120041: iFruchter et ail 120061: 
iSavaglio et al.ll2009HSvensson et al.ll201Q[ ). and late time 
bumps rese mbled supernova componen ts in light curves 
of LGRBs (|Zeh et al.ll200l 120051 l2006h . 

The observational association between GRBs and su- 
pernovae has provided strong support to a scenario, so- 
called collapsar model, in which LGRBs are assumed to 
be originated i n the collapse of a massive stellar core 
to a b lack hole ()Wooslevlll993l ). iMacFadven fc Wooslevl 
(|1999D outlined possible scenarios of driving LGRBs. 
In the collapsar model, a central core of a mas- 
sive star is required to be rotating rapidly enough 
that a massive accretion disk can be formed around 
a black hole. Then, pair annihilation of neutrinos 
emitted from the accretion disk to electron-positron 
pairs could s upply sufficient en e rgy to induce relativis- 
tic outflows (lEichler et al Ill989t iMeszaros fc Reeslll99^: 
iNaravan et al.lll992t iMochkovitch et al.lll993l ). The rel- 
ativistic outflows are expected to form a GRB fireball. 
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In addition, it is suggested that strong magnetic fields of 
order 10 15 G, if they are present, coul d play an active 
role in driving the relativistic outflows flN akamur a et al.l 
[19921 iNarayan et alJll99l iLvuikov ll2006l ). 

There are three possible varie ties in collapsar model 
(iHeger et al.l 120031 ): In Type I (iMacFadven fe Wooslevl 
11999J) and Type II (jMacFadven et al.l 1200 ID collapsar 
models, a proto-neutron star is assumed to be formed 
initially and a shock wave is launched. Then, in the Type 
I collapsar, the proto-neutron star collapses promptly 
to a black hole because the shock wave is weak, while 
in the Type II collapsar, a black hole is formed by a 
fallback process long after the proto-n eutron star forma- 
tion. In the Typ e III collapsar model ()Heger et al.ll2003l : 
iFrver et al.ll200l[ ) , a black hole is directly formed without 
formation of proto-neutron star. 

Recently, two LGRBs (GRB060505 and GRB060614) 
which are not li kely to be accompanied by a supernova 
were discovered dFvnbo et al.ll2006t iGehrels et al.l 120061: 
iGal-Yam et al]|2006l: iDella Valle et al.ll2006al) . The host 
galaxy of GRB0600505 is a star-forming galaxy similar 
to that of canonical LGRBs. Such LGRBs might be as- 
sociated with the Type I or Type III collapsar. Note 
that there is d ebate about the lack of supernova fe ature 
in GRB06014 (ICobb et al.1 [20061: IDado et al.l[2008h and 
it has been discussed that the duration of G RB060505 is 
abou t 4 second and it may be a short GRB (Ofek et al. 
[2001 . 

Because the observed supernovae associated with 
LGRBs are Type Ib/c an d the relativistic jets ha ve to 
reach the stellar surface (jZhang fc Wooslevl I2004D . the 
progenitors should have lost their envelope before the on- 
set of stellar core collapse; otherwise a peculiar evolution 
path is required. Due to these reasons, the progenitors of 
LGRBs are now believed to be rapidly rotating massive 
Wolf-Rayet (WR) stars. However, ordinary WR stars are 
known to be accompanied by strong stellar winds driven 
by radiation pressure which lead to a rapid spin-down of 
the stellar core. Here, a serious problem concerning col- 
lapsar model is that according to stellar evolution calcu- 
lations, it is very difficult to produce pre-collapse cores 
which satisfy both the requirement of collapsar model 
and the association of Type Ib/c supernova, if magnetic 
torque s and standard mass-los s rates are taken into ac- 
count (|Wooslev fc Heger|[2006l) . 

To resolve the ab ove dilemma , seve ral models have 
been proposed (se e IFrver et all (120071) for a review). 
Ilzzard et al] ([20041 ) and iPodsiadlowski et al] (j2004f ) pro- 
posed binary-interaction models, in which the tidal force 
in a close binar y keeps a helium star i n sync hronous, 
rapid rotation. Ivan den Heuvel fc Yoonl ([20071 ) showed 
that a helium star in a close binary with a compact com- 
panion (i.e., neutron star or black hole) can retain suffi- 
cient angular moment um to form a progenitor of a GRB. 
IFrver fc Hegerl ([20051 ) suggested a binary-merger model 
and showed that a merger of two helium cores during the 
common-envelope inspiral phase can produce a rapidly 
rotating core which satisfies the requirement of the col- 
lapsar models. 

On the o t her hand, Yoon fc Langer (jYoon fc Langerl 
l200l [20061 lYoon et al.1 [2001) and Woosley & Heger 
I Wooslev fc Hegerl 120061 ) recently showed that a single 
star can fulfill the requirements of the collapsar models 
if it is initially rapidly rotating (> 50% of the Keplerian 



velocity at the equatorial surface) and of low metallic- 
ity (Z/Zq < 0.1). Note that the low metallicity could 
keep the stellar radius small er and also reduce the mass 
loss (j Wooslev fc Hegerl 120061 ). Both effects suppress the 
loss of angular momentum from the star. The rapid ro- 
tation results in a short mixing timescale, which could 
help achieving a chemically homogeneous state through- 
out the hydrogen burning phase. In this case, a single 
star could become a rapidly rotating WR star without 
losing the hydrogen envelope through the stellar wind, 
avoiding the red giant phase that otherwise would cause a 
significant decreas e of the core angular m omentum due to 
magnetic torques (jYoon fc L anger 2006|) . It is also noted 
that the chemically homogeneous evolution is likely to 
occur for the tidally spun-up star in a binary system 
(jGantiello et al.ll27)07l ). 

There are several supports to the chemically- 
homogeneous-evolution model. Recent observations have 
indicated tha t LGRBs may prefer a low metallicity 
environment dFruchter et al.l [2001 iStanek et~aT1 120061: 
iModiaz et all 120081: ISvensson et al.ll2010D . If the binary 
merger model resulted in most of the LGRB progenitors, 
such dependency would not be found. 

Gravitational collapse of population III (Pop III) stars, 
which are assumed to be formed from metal-free gas, 
may be accompanied by LGRB at a very high red- 
shift (|Schneider et al.l [20021: iBromm fc l.oo"fl [21 )()<>:■. Nu- 
merical simulations have suggested that Pop III stars 
would be predom i nantly very massive with M > lOOMm 
(lOmukai fc Para! I200ll 120031: iNakamura fc Umemural 



2001: lAbeT et al.l 2002: Br omm et al. 2002). Such a mas- 
sive star may collapse directly to a black hole without 
producing supernova explosion (Type III collapsar). 

In addition, an attempt to constrain th e characteristics 
of LGRB progenitors has been made by iCampana et alJ 
(2008), who studied in detail an absorption pattern in 
the X-ray spectrum of GRB060218 and found an ex- 
tremely low O /N ratio in the surrounding of the progen- 
itor, reaching a conclusion that only a progenitor star 
characterized by a fast rotation and subsolar metallicity 
could explain it. 

All of the above progenitor models of LGRBs are 
anomalous in the sense that they arc different from the 
progenitors of ordinary supernovae. Qualitatively speak- 
ing, the progenitor models should produce a core of larger 
angular momentum than the ordinary supernova cores. 
Also, the central entropy of the core would be higher than 
the ordinary supernova cores because of its high mass: 
The chemically homogeneous models tend to predict a 
well-mixed, larger core with higher central entropy than 
the ordinary supernova core. It is also expected that the 
object formed after the binary merger will have a higher 
entrop y, if the mass ratio o f merging stars is not far from 
unity (jSuzuki et afl 120071 : iGaburov et"aTI 120081 ) . Thus, 
LGRB progenitor cores may be modeled by a rapidly 
rotating, higher-entropy core, regardless of their forma- 
tion processes. Based on this assumption, in this paper, 
we perform collapse simulations of a very massive stellar 
core with a fairly high value of entropy (s = 8ks per 
baryon) to study effects of higher-entropy. 

A number of hydrodynamic simulations have been 
performed for studying gravitational collapse of such 
rapidly rotating, higher-entropy core in the context of 
collapsar model: for the Type I collapsar model, see 



3 



iProga et all (12003ft .iFuiimoto et al J (12006ft. iDessart et all 



(200S 



Nakataki 



lLopez-Camara et alJ 



(2009), 
(2009), 



lHarikae et al 
and TotT et alJ 



for th e Type II collapsar model, see iMacFadve n et al 



(2009) 



(2001); for the Type III col l apsar model, see iFrver et al 
( 2001 ) IShibata fc Shapirol (l2002h. ISekiguchi fc Shibatal 
( 2007) , ISuwa et alJ (I2007bft . and lLiu et al.l (120071) ). Most 
of the simulations were perfo rmed in the New t onian 
or pseudo-Newtonian gravity ( MacFad ven et al.| 



200C 



Frver et al|2001l:lProga et all2003MFuiimoto et al.l 
Suwa et all 12007b!: iDessart et alJ 120081 : lHarikae et alJ 



2001; 



20091: lLopez-Camara et £01 12009) . In such simulations, 
inner regions of core (r < 5-20rs where r§ is the 
Schwarzschild radius) are excised, and consequently, 
increase of the overall efficiency of accretion according 
to the black hole spin from « 6% (zero spin) to « 42% 
(maximal spin) cannot be taken into account. The 
black hole spin has significant effects on structure of 
the accretion disk, because it dramatically changes the 
spacetime metric near the black hole, where most o f 
accretion power is released (|Chen fc Beloborodovl l2007). 

Also, to guarantee formation of a centrifugally sup- 
ported accretion disk at radii larger than the excised 
radius, most of Newtonian studies adopted angular mo- 
mentum distributions that are well above the threshold of 
the disk formation: The specific angular momentum j for 
a large fraction of the core is assumed to be much larger 
than that at the innermost stable circular orbit (ISCO), 
iisco- In such cases, gravitational energy will not be ef- 
fectively converted into thermal energy due to the large 
radii. Rather these models rely on subsequent hypotheti- 
cal viscous h eating for generating large am ount of energy. 
By constant. iLee fc Ramirez- Ruiz! ([20 06) performed sim- 
ulations of low angular momentum accretion flows into 
a black hole in the Newtonian framework. They found 
that a thin accretion disk is formed for j < 1.9rgc 
while a thick torus i s form ed for j < 2.1r s c (see also 
lLopez-Camara et alJ (|2009ft ). lHarikae et all ([2009ft also 
found similar results. 

To self-consistently follow formation of a black hole and 
a surrounding disk, a fully general relativistic simulation 
for the collap se of rapidly rotating mass ive star was first 
performed bv lShibata fc Shapirol ([2002ft . Unfortunately, 
they could not follow the subseque nt evolution of an ac- 
cretio n dis k around the bla ck hole. ISekiguchi fc Shibatal 
(|2007ft and lLiu et all ()2007l ) performed fully general rel- 
ativistic simulations of collapsar, successfully following 
formation of an acc retion disk and a n early evolution of 
the disk. Recently, lOtt et all ([2011ft performed simula- 
tions in the context of the collapsar scenario an d extract- 
ing th e gravitational wave signature from it. iNakatakl 
(2009) performed a long-term general relativistic simula- 
tion in a fixed Kerr black hole background. However, in 
these general relativistic simulations, relevant microphys- 
ical processes such as neutrino cooling were not taken 
into account. 

In this paper, we for the first time report the re- 
sults of fully general relativistic simulations for the col- 
lapse of a rapidly rotating, high-entropy core, taking 
into account detailed microphysics; a nuclear-theory- 
based finite-temperature equation of state (EOS), weak 
interaction processes such as electron capture and pair- 
neutrino processes, and neutrino cooling. We focus on 
self-consistently clarifying the formation process of a ro- 



tating black hole and surrounding accretion disk, and 
subsequent long-term evolution of this system. We will 
show how the black hole is formed and evolved, and also 
clarify the physical condition for the disk or torus in the 
vicinity of the black hole. In particular, this is the first 
work that clarifies the geometrical structure, thermal 
property (such as chemical composition, chemical poten- 
tials, and entropy), neutrino optical depth, and neutrino 
luminosities of the accretion disk in the framework of full 
general relativity. 

The paper is organized as follows. We first briefly 
summarize the basic equations, the input physics, and 
numerical setup in Section [21 The main results are de- 
scribed in Section [3] Discussion of our results together 
with prospects for GRB production are given in Section 
[31 Section [5] is devoted to a summary. Throughout this 
paper, h, fcs, c, and G denote the Planck's constant, 
the Boltzmann's constant, the velocity of light, and the 
gravitational constant, respectively . W e a dopt the geo- 
metrical unit c = G = 1 in Sections 12.11 and 12.21 which is 
commonly used in numerical relativity. 

2. SETTING 

2.1. Einstein's equation and gauge conditions 

The standard variables in the 3+1 decomposition of 
Einstein's equation are the three-dimensional metric 7^ 
and the extrinsic curvatu re Kg on th e three-dimensional 
hypersurface denned by ([York l Fl979) 



7*w = 



- g^v 1 fi^nis . 



where g^ w is the spacetime metric, 



(1) 
(2) 

is the unit normal 



to a three-dimensional hypersurface, and £, n is the Lie 
derivative with respect to the unit normal n M . Then we 
can write the line element in the form 



ds 2 = 



2 alt 2 + j ij (dx i + P l dt){dx 3 + (3 j dt), 



(3) 



where a and are the lapse function and the shift vector 
which describe the gauge degree of freedom. 

Numerical simulatio n is performed in the 
BSSN formulation (IShibata fc Nakamural 119951: 
iBaumgarte fc Shapiro! 11999ft hi which the spatial 
metric 7y is conformally decomposed as 7y = e 4 ^jij 
where the condition det(7y) = 1 is imposed for the 
conformal metric jij . From this condition, the conformal 
factor is written as <j> — j^ln7 and 7 = det(7ij). The 
extrinsic curvature Kij is decomposed into the trace part 
K and the traceless part A^ as = A. L j + (l/2>)^ijK 
The traceless part Ay is conformally decomposed 
as A^ = e^Aij. Consequently, the fundamental 
quantities for the evolution equation are now split into 

7y, K, and Ay. Furthermore, the auxiliary variable 
Fj = S^dkla is introduce d in the BSSN formulation 
(|Shibata fc Nakamuralll995ft . 

To stably follow the spacetime after appearance of a 
bla ck hole, we evolve W = e~ 2 ^ instead of 4> follow- 
ing iMarronetti et al.1 (|2008ft . The primary reason is that 
<j> diverges at the center of a black hole in the vertex- 
center grid. With the choice of W, such pathology can b e 
avoided, as first pointed out bv lCampanelli et all (|2006ft . 
in which \ = e _4< ^ was used instead of W. Merits of 
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using W are that (i) the equation for the Ricci tensor is 
slightly simplified, (ii) no singular term appears in the 
evolution equations even for W -> 0, and (iii) the deter- 
minant of 7ij is alway s positive ( Marron etti et all 120081 : 
lYamamoto et al.ll2008f) . 

We assume axial and equatorial symmetrie s of the 
spacetime and the so-called Carto on method (jShibatal 
120001 12003at lAlcubierre et al.ll200lD is adopted to avoid 
possible problems around the coordinate singularities of 
the cylindrical coordinates. In the present code, we use a 
4th-order finite difference scheme in the spatial direction 
and a 3rd-order Runge-Kutta scheme in the time integra- 
tion. The advection terms suc h as P l dj(l> are evaluat ed 
by a 4th-order upwind scheme (Briigman n et al.1 l2008) . 

As the gauge conditions for the lapse, we us e a dynam- 
ical slicing f cf. lAlcubierre fc' Briigmann 2001): 

d t a = -2Ka. (4) 

It is known that this dynamical slicing enables to perform 
a long-term evolution of neutron stars as well as has a 
strong singularity avoidance property in the black hole 
spacetime. The shift vector is determined by solving the 
following dynamical equation (Shibata 2003b) 



d t /3 k = ^ kl (Fi+Atd t F l ). 



(5) 



Here the second term in the right-hand side is neces- 
sary for numerical stability, and At denotes the numeri- 
cal timestep. 

2.2. Hydrodynamic equations coupled to general 
relativistic leakage scheme 

Recently, Sckiguchj (|2010al lH) developed a fully general 
relativistic hydrodynamic code implementing a nuclear- 
theory-based finite-temperature EOS, self-consistent 
electron and positron captures, and neutrino cooling by 
a general relativistic leakage scheme. Neutrino heating 
is not included in the current version of leakage scheme. 
Since we assume the axial and equatorial symmetry of the 
spacetime, the hydrodynamics equations are solved in the 
cylindrica l coordinates (w, (p, z) where w — \J x 2 + y 2 . 
We follow Sckiguchi (2010b) for a solution of the hydro- 
dynamic equations to which the readers may refer for 
the details. In this section, we adopt the geometrical 
unit c = G = 1. 

2.2.1. Energy-momentum conservation equation 

The basic equations of general relativistic hydrody- 
namics with neutrinos are 

Vo(T TbtaIja = Vq [ (T Fja + (j-)^ = Qj (6) 

where (T Total ) Q ^ is the total energy-momentum ten- 
sor, and (T F ) a p and (T v ) a p are the energy-momentum 
tensor of f luids a nd neutrinos, respectively. Following 
Sckiguchi (2010b), the neutrino energy- momentum ten- 
sor is decomposed into 'trapped-neutrino' ((T"' T ) a p) and 
'streaming-neutrino' ((T' l/ ' S ) Q/ 3) parts as 



u.S\ 



(7) 



Here, the trapped-neutrino part phenomenologically rep- 
resents neutrinos which interact sufficiently frequently 
with matter, and the streaming-neutrino part describes a 



phenom enological flow of neutrino s streaming out of the 
system. iLiebendorfer et all (|2009T > developed a more so- 
phisticate method in terms of the distribution functions 
of trapped and streaming neutrinos in the Newtonian 
framework. 

Streaming-neutrinos are produced with a leakage rate 
Q^ ak , according to 



v (r s ) a 







-)lcak 



(8) 



On the other hand, the trapped-neutrino part is com- 
bined with the fluid part as 

T a p = {T F ) aP + (T"> T ) Q/3 . (9) 

Then the equation for T a p is 

V a T^ = -Q^ ak . (10) 

We solve Eqs. ([5]) and (fit)]) for the energy-momentum 
conservation equation. 

The energy-momentum tensor of the fluid and trapped- 
neutrino parts {T a p) is treated as that of the perfect fluid, 



T a p = {p + pe + P)u a Uf) + Pg a/ : 



(11) 



where p and u a are the rest mass density and the 4- 
velocity. The specific internal energy density (e) and 
the pressure (P) are the sum of contributions from the 
baryons (free protons, free neutrons, a-particles, and 
heavy nuclei), leptons (electrons, positrons, and trapped- 
neutrinos), and photons as, 

(12) 
(13) 



P = P B +Pe+Pu+Pph, 

£ = £b + £e + £i> + £ph , 



where subscripts 'B\ 'e', 'pft.', and V denote the compo- 
nents of baryons, electrons and positrons, photons, and 
trapped-neutrinos, respectively. 

The streaming-neutrino part, on the other hand, is set 
to be a general form of 



{T v - S ) a p = En a np + F a np + F p n a + P a/3l 



(14) 



where F a n a — P a pn a — 0. In order to close the system, 
we need an explicit expression of P a /3- In this paper, we 
adopt a simple form P a p = xEj a /3 with x = 1/3- Then 
we solv e Eq. (|8"|) in a h igh resolution shock capturing 
scheme (|Sekiguchill2010bD . 

The closure relation employed in this paper is not very 
physical. Also, recall that we do not consider the so- 
called neutrino heating in this paper. To treat the neu- 
trino heating accurately, a more sophisticated closure re- 
lation is required. However, such a study is beyond the 
scope of this paper. A more sophisticated treatment of 
neutrino transport equations, together with incorporat- 
ing t he neutrino heating , will be needed in the future 
(e.g.. lShibata et al.l[20rl . 

2.2.2. Lepton-number conservation equations 

The conservation equations of the lepton fractions are 
written schematically as 

dY e 
dt 

dY Uc 



"7e 



(15) 
(16) 
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dt 
dY Vx 
dt 



(17) 
(18) 



where Y e , Y Ue , Yp e , and Y Vx denote the fractions per 
baryon number for electrons, electron neutrinos, elec- 
tron anti-neutrinos, and p and r neutrinos and anti- 
neutrinos, respectively. Here we consider, as local 
reactions, the electron capture, the positron capture, 
electron-positron pair annihilation, plasmon decay, and 
the Bremsstrahlung radiation of pair neutrinos, where 
v and v denote the three flavors of neutrinos and anti- 
neutrinos. 
The source terms are given by 



7e 




_ local 


(19) 


ll>e 


— ^local 

IU e 


_ ^l cak 
iv c ' 


(20) 


7p s 


local 

= 7 Pe 


leak 

~7p e , 


(21) 


1v* 


— local 




(22) 



where <-y local 's and <y ieak 's are the local production and 
leakage rates of each species of neutrinos, respectively. 
Because <y^ ocal >g are characterized by the timescale of 
weak-interaction processes t wp ~ \Y e /Y e \ which can 
be much sho rter than the dynamical timescale (e.g., 
lBruennll"l985D . a straightforward explicit solution of Eqs. 
(fT5*]) - ([T£t leads, in general, to a numerical i nstability. 
Therefo re we follow the procedure proposed in lSekiguchil 
(|2010b| ) to solve the equations stably in an explicit man- 
ner. 

First, in each timestep n, the conservation equation of 
the total lepton fraction (Yj =Y e — Y„ e + Yp e ) , 



dt 



(23) 



is solved together with the conservation equation of Y v% , 
Eq. (fT8"|) . in advance of solving the whole of the lepton 
conservation equations (Eqs. (TTBl - (HBJ). Then, as- 
suming that the /3-equilibrium is achieved, values of the 
lepton fractions in the /3-equilibrium (Y/, Yj? , and Y^) 
are calculated from the evolved value of Y . 

Second, regarding Y& and Y$ as the maximum allowed 
values of the neutrino fractions in the next timestep n+1, 
the source terms are limited so that each value of Y„'s 
in the timestep n + 1 cannot exceed that of Y$ 's. This 
limiter procedure enables to solve explicitly the whole of 
the lepton conservation equations (Eqs. (fT5|) - (fT8|l). 

Third, the following conditions are checked, 



u p + n e < fi n + fj, Ve . 
[J, n - n e < p p + p, De 



(24) 
(25) 



where p, p , p n , p e , p, Ue , and p„ e are the chemical potentials 
of protons, neutrons, electrons, electron neutrinos, and 
electron anti-neutrinos, respectively. If both conditions 
are satisfied, the values of the lepton fractions in the 
timestep n + 1 are set to be those in the /3-equilibrium 
value; Yf , Yjf , and Y^ . On the other hand, if either or 
both conditions are not satisfied, the lepton fractions in 
the timestep n + 1 is set to be those obtained by solving 
the whole of the lepton-number conservation equations. 



2.3. Microphysics 
2.3.1. Equation of state 

I n this paper we employ a tabulated EOS derived 
by IShen et all (|1998f ) , which is based on the Briickner- 
Hartree-Fock-type relativistic mean field theory. The 
maximum gravitational mass of a cold spherical neutron 
star in this EOS is much larger than the canonical neu - 
tron star mass 1AM Q as w 2.2M (jShen et al.lll998[) . 
The framework of the relativistic mean field theory is ex- 
tended with the Thomas-Fermi spherical cell model ap- 
proximation to describe not only the homogeneous mat- 
ter but also an inhomogeneous one. 

The thermodynamical quantities of dense matter at 
various sets of (p,Y p ,T) are calculated to construct 
the numerical data table for simulation. Here Y p is 
the total proton fraction per baryon number. The 
original table covers a range of density 10 51 -10 15 4 
g/cm 3 , proton fraction 0.0-0.56, and temperature 0— 
100 MeV, which are required for supernova simulation. 
The original table has been e xtended to higher density 
(jSumivoshi et al.l 120071 I2008D and higher temperature 
(jNakazato et all I2008D ranges of 10 51 -10 17 g/cm 3 and 
0-400 Me V, which are required f or following black hole 
formation (jSumivoshi et al.|[2006l) . 

It should be noted that the causality is guaranteed 
to be satisfied in this framework, whereas the sound 
velocity sometimes exceeds the velocity of the light 
in the non-relativistic fra mework, e.g., in the EOS by 
lLattimer fc Swestvl (|1991[ ) . This is one of the benefits of 
the relativistic EOS. 

To consistently calculate the pressure and the internal 
energy of electrons and positrons, the charge neutrality 
condition Y p = Y e should be solved to determine the elec- 
tron chemical potential \x e for each value of the baryon 
rest mass density p and the temperature T in the EOS 
table. Namely, it is required to solve the equation 



( r r\ - P Ye 
n e (fi e , 1 ) = n_ - n + = 

m„ 



(26) 



in terms of \i e for given values of p, T, and Y e {= Y p ). 
Here, m u = 931.49432 MeV is the atomic mass unit, and 
n_ and n + are the total number densities (i.e., including 
electron-positron pairs) of electrons and positrons, re- 
spectively. Then, assuming that electrons and positrons 
obey the Fermi-Dirac distribution, the number density, 
the pressure, and the internal energy density of electrons 
and positrons are calculated in a standard manner (e.g., 
ICox fc GiulilH968l) . 

The pressure and the specific internal energy density 
of photons are given by 



a r T 



P 



(27) 



where a r = (tt 2 kg) / (15c 3 h 3 ) is the radiation constant. 

In this paper, trapped-neutrinos are assumed to in- 
teract sufficiently frequently with matter that be ther- 
malizcd. Therefore they are described as ideal Fermi 
gases with the matter temperature. From the numeri- 
cally evolved neutrino fractions YJ 5 ™ 1 , the chemical po- 
tentials of neutrinos (p u ) are calculated by solving 



Y;™ 1 = Y^ V .T) = ^n v {p Vl T). 

P 



(28) 
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Then the pressure and the internal energy of trapped- 
neutrinos are calculated in the same manner as for elec- 
trons, using fi u and the matter temperature. 

2.3.2. Weak interaction and leakage rate 

Following ISekiguchil (|2010b| ). the leakage rates are de- 
fined by 



(1- 



diff 



leak 



-br 



local 



+ e 



e bTv Q* 

-67V local 

IV i 



(29) 
(30) 



where t v is the optical depth of neutrinos and b is a pa- 
rameter which is typically set as b^ 1 — 2/3. The optical 
depth can be computed from the c ross sections follow- 
ing an often employed prescri ption ([Ruffert et al.l[l996l : 
iRosswog &: Liebend orfcr 2003|) : The optical depth is cal- 
culated by 

t v = min[r^ ,t*,t£] , (31) 

where t™ ', r*, and r£ are the optical depths along zu, z, 
and the radial directions, respectively. We calculate, for 
example, t£ by 



T*(07,JS) = 



k v {w 1 z')dz' , 



(32) 



where n v is the opacity and z out denotes the outer bound- 
ary in the z-direction. and r£ are calculated in a 
similar manner. 

Then, because Q l ° ak should be regarded as the emis- 
sivity of neut rinos measured in the fluid rest frame. P 1 ' ilk 
is defined as (jShibata et al.ll2007t ISekiguchil l2"010al|t 



^lcak 



^lcak^ 



(33) 



As the local production reactions of neutrinos, we con- 
sider the electron and positron captures (7®° and 7p e c ) fol- 
lowing ^lliZiES] (|1985j ). the electron-positron pair an- 
nihilation ('Jue^e ^ or electr on- type neutrinos and 7^ "„ for 
the other type) following iCooperstein et ail (|1986t ). the 
plasmon decays (7^^ and 7£Jp*,) following iRuffert et al.l 
(1996), and the Bremsstrahlung processes (7^ e ms and 

7^r S ) following iBurrows et al.l ()2006[ ). Then, 'the local 
reaction rates for the neutrino fractions are 



local = oc , pair plas Brcms 

/ I/ e / U a { I h> a V e 1 / Z-'e Z-'e 1 / U e V e ' 

local _ pc , pak plas Brems 



local _ pair 



plas 



Brcms 



(34) 
(35) 
(36) 



Similarly, the local neutrino energy emission rate Q l ° cal 
is given by 



glocal = Q cc + gpc + 3 (Q pan + g pras + 



-jplas 



Qplas 



jBrcmsj 
nBrcms^ 



(37) 

(da 



The explicit forms of t he loca l rates in Eqs. (f34|) 

are found in iSekiguchil (|2010bl ) 

We follow the recent work bv IRosswog fc Lie bendorfcr 
(2003) for the diffusive neutrino emission rates 7^ lff and 
Qf s in Eqs. ([29]) an d fl5DJ|. The explic it forms of 7^ iff 
and Qf B are found in ISekiguchil (|2010bft . 
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Fig. 1. — Radial profiles of density (upper panel) and temperature 
(lower panel) of the initial configuration. 
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Fig. 2. — Radial profiles of specific angular momentum for the 
slowly, moderately, and rapidly rotating models. 



2.4. Initial model 

Because there are no realistic models of rotating pro- 
genitors derived by multi-dimensional pre-collapse evo- 
lution calculations or no binary progenitor models, we 
prep are approximate initia l models in the following man- 
ner (|Nakazato et al.l |2007f) . We first calculate a spher- 
ical equilibrium configuration with a constant electron 
fraction of Y e — 0.5 and with a constant entropy per 
baryon s = 8fcs. We set the central density to be 
p c w 10 s g/cm 3 . The corresponding central tempera- 
ture is T c « 9 x 10 9 K, which is higher than the critical 
temperature for t he photo-dissociation of heavy nuclei to 
occur. Following Nak azato et al.l (j2007l ). we define the 
outer boundary of the 'iron core' to be where the tem- 
perature is 5 x 10 9 K. Note that most of heavy nuclei in 
inner parts of this 'iron core' in fact are already photo- 
dissociated. Then the mass and the radius of the core are 
A/iron ~ 13M Q and n ron 7000 km. In numerical simu- 
lation we follow a region of r tot m 14000 km (> rj ron ) in 
which the total mass of M tot « 23M© is enclosed. The 
radial profiles of density and temperature are shown in 
Figure [TJ 

For the purpose of reference, we note that our initial 
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Fig. 3. — Distributions of the averaged specific angular momen- 
tum for slowly (blue curve), moderately (green curve), and rapidly 
(red curve) rotating models. The specific angular momentum re- 
quired to support a fluid element in a circular orbit at ISCO around 
a Schwarzschild black hole and a maximally rotating Kerr black 
hole of mass m»(j) is shown together (black dotted curves). The 
blue triangles, green squares, and red circles indicate the numerical 
results for the path along which specific angular momentum and 
mass of the black hole formed in the collapse of the slowly, mod- 
erately, and rapidly rotating models follow, respectively (see Sec. 
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Fig. 4. — The spin parameter distribution, q*(j), and the specific 
angular momentum at ISCO, Jisco(j)/-^*> m units of c/G, where 
M* is the total baryon mass. 

model might correspond to entropy p er baryon for a sta r 
with initial mass of » 1 20-130Mp) dBond et all 11984ft . 
However, a recent study ()Waldm an 2008) predicts that 
such massive stars will undergo a pulsational pair insta- 
bility and considerable mass loss, resulting in hydrostatic 
degenerate iron cores of mass with ~ 3Af Q , which is dif- 
ferent from the initial mod el adopted in this p aper. The 
300 Mq progenitor used bv lFrver et al.l (|2001ft has a cen- 
tral entropy of ~ 8fcs per baryon. However such a very 
massive model does not form an iron core in hydrostatic 
fashion, but rather goes unstable much earlier burning 
phase. Note that these are results in a spherical sin- 
gle star with solar metallicity. Anomalous stars, such as 
stars in interacting binar y and Pop III stars, m ight form 
such high-entropy cores ( Nakazato et al.|[2007ft . 
Little is also known about the angular momentum dis- 



tt(w) = f2fj exp 



1 



2 {w 2 + Rl) 



Thus, we employ the 



(38) 





" w 2 ' 


exp 


R 2 



where w — y 1 x 2 + y 2 , and ft, Ro, and R c are parameters 
which control the degree of differential rotation. The ex- 
ponential cut-off factor is introduced by a practical rea- 
son for numerical simulation: if the specific angular mo- 
mentum in the outer region of the core is too large, the 
matter escapes from the computational domain. How- 
ever, the most part of the 'iron core' is almost uniformly 
rotating. We fix the values of Ro and R c as Ro = r tot /5 
and R c = r tot /8, respectively. We vary fi as 0, 0.4, 0.5 
and 0.6 rad/s (hereafter referred to as spherical, slowly 
rotating, moderately rotating, and rapidly rotating mod- 
els). The rotation period in the central region is « 10- 
15 s. This is by one order of magnitude longer than the 
dynamical timescale (Gp c ) -1 / 2 ~ 0.4 s. Thus, the pro- 
genitor star is not assumed to be rapidly rotating. The 
profiles of specific angular momentum along the cylindri- 
cal radius are plotted in Figure [2J 

Figure [3] plots an averaged specific angular momen- 
tum distribution defined by J*(j)/m*(j). Here, j is the 
specific angular momentum of a fluid element, which is a 
conserved quantity in axially symmetric spacetime in the 
absence of viscosity. m*{j) is a rest mass distribution as 
a function of j , which is the integrated baryon rest mass 
of fluid elements with the specific angular mom entum 
less than j, defined by (IShibata fc Shapiroll2002ft 

m*(j) = 2tt / p t r 2 drd{cos9). (39) 
Jj'<i 

Similarly, J*(j) is an angular momentum distribution de- 
fined by 



p*j'r 2 drd(cos 9). 



(40) 



<3 



These conserved quantities are often used in general 
relativistic st udy to predict a possi b le outcome of 
the collapse (IShibata fe Shapiro! 120021: iShapiro I l2004t 
ISekiguchi fc Shibatall2004ft~ 

It should be noted that the specific angular momen- 
tum considered in this paper is rather small for a large 
fraction of fluid elements, in the sense that it is smaller 
than the angular momentum required for a fluid ele- 
ment to stay outside the innermost stable circular or- 
bit (ISCO), jisc0 7 around a Schwarzschild black hole. 
In this sense, our model is 'sub-Keplerian'. This is by 
contrast with many of previous models in which the 
specific angular mo mentum of well above 7 t sc:o i s usu- 
ally imposed (e.g.. iMacFadven fc Wooslevl (119991), but 
see | Lee fc Ramirez-Ruiz! (12006ft . iLopez-Camara et al.l 
(2009). and lHarikae et all (|2009ft ). In the present con- 
dition, the fluid elements of such small specific angular 
momentum form a black hole, while those of large specific 
angular momentum does a disk (torus). 

Now, to infer the evolution of a black hole surrounded 
by accreting materials, let us consider ISCO around a 
hypothetical black hole located at the center. If the value 
of j of a fluid element is smaller than that at the ISCO, 
jisco , f°r the hypothetically formed black hole, the fluid 
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element will fall into the seed black hole eventually. The 
value of jisco will change as the ambient fluid elements 
accrete into the black hole. If jisco increases as a result 
of the accretion, more ambient fluid elements will fall 
into the black hole. On the other hand, if jisco decreases 
during the accretion, the accretion into the black hole will 
be suppressed, and then, the black hole will approach to 
a quasi-stationary state with a small accretion rate. 

To estimate the value of jisco > we assume that the 
spacetime metric can be instantaneously approximated 
by that of a Kerr spacetime of mass m*(j) and the non- 
dimensional spin parameter q*(j) = cJ*(j)/Gm*(j) 2 . 
On these appro ximations, we may comput e jisco of a 
black hole (e.g.. iShapiro fc Teukolskvlll983f ). 

For all the models considered in this paper, g*(j) is 
smaller than unity for a fraction of fluid elements with 
small specific angular momentum. As a result of this fact, 
these fluid elements can form a black hole in the dynam- 
ical timescale. However, this will not be the case for the 
initial condition with g*(j) > 1 in an inner region. In this 
case, a black hole will not be formed directly because the 
Kerr space time with the spin parameter greater than 
unity contains a naked singularity . Instead, a rotating 
oblate object will be the o utcome ((Saiio fc Hawkel [20091 : 
iSekiguchi fc Shibatal I2004H . Such an oblate object will 
be unstable against nonaxisymmetric deformation, and 
then, angular momentum will be transported by the hy- 
drodynamic torque from the inner region to the outer 
one. As a result of a sufficient amount of angular mo 
mentum transpor t, a black hole will be eventually formed 
(jZink et al.1 120071 ). This suggests that the timescale for 
black hole formation may be determined by the timescale 
for the angular momentum transport. We do not con- 
sider this possibility in this paper. 

Figure 2] plots the spin parameter distribution (<7*(j)) 
and jisco (j) = jiSCO [m*(j),q*(j)] as functions of m*(j). 
This figure clearly indicates that the value of jisco (j) 
takes the maximum at ™*(j) ~ 12,16, and 20Mq for 
the rapidly, moderately, and slowly rotating models, re- 
spectively. These values show a possible final value of 
black hole mass, which is smaller than the total mass of 
the system. This indicates that a certain fraction of the 
material with mass > M© will form a disk around the 
black hole. It should be noted that the curves of Fig- 
ures |3] and E] indicate the possible evolution path of the 
black hole only approximately. In determining jisco as 
a function of m*(j), we assume that a fluid element of 
smaller value of j falls into black hole earlier. However, 
this is not always the case in the dynamical evolution of 
the system, because the material in the outer region near 
the rotation axis has a small value of j and falls into the 
black hole in a late time. 

2.5. Analysis of black hole and accretion disk 

The formation of a black hole is ascertained by finding 
apparent horizon (Shibata 1997). Then, we calculate two 
geometrical quantities which possibly characterize mass 
of a black hole. One is an irreducible mass defined by 



M irr = 



.4 



H 



G V 16tt' 



(41) 



where Ah is the area of the apparent horizon. The other 
mass is associated with the circumference proper length 



along the equatorial surface C e : 

M r 



2 C e 



c 



(42) 



This should agree with the mass of a Kerr black hole in 
the stationary axisymmetric spacetime. Note that in the 
case of a Schwarzschild black hole M- m = M ce . 

We also estimate black hole mass using an approximate 
conservation law, 

M con — A^ADM — M*,r>r AH i (43) 

where Madm is the ADM mass of the system and 
M*, r>rAH is the rest mass of baryons located outside the 
apparent horizon. It is suggested that M ce may be a 
good indicator of mass of a black hole even in the pres- 
ence of a massive accretion disk ([Shibatal 120071 ). As we 
shall see in Section |3l M ce and M con agree approximately 
with each other, and thus, we use M ce as the black hole 
mass, namely, 

M BH = Mce W Mcon- (44) 

The non-dimensional spin parameter q of a Kerr black 
hole can be calculated from the ratio between polar and 
equatorial circumferential radii of event horizon, C p and 
G e , 




where f + = 1 + \f\ 
Kerr black hole, 



(45) 



The definition of M\ rr for a 



Mi 



BH 



(46) 



may be also used to estimate the black hole spin. How- 
ever, by contrast with M ce , C p /C e and M; rr /MBH are 
not very good indicator s of the black hole spin when 
a massive disk presents (jShibatal 120071 ) . In the case of 
equilibrium configuration of a black hole surrounded by 
a massive disk, it was found that a spin parameter esti- 
mated by Eqs. (|45|) and (|46|) decreases with the increase 
of disk mass and with the decrease of the inner edge of 
a disk. Accordingly, a spin parameter estimated by Eqs. 
(j45|) and (|46|) may contain an error of Ag ~ 0.1 because a 
massive accretion disk falling into a black hole is formed 
in the present study. 

We note that we approximately calculate C p , C e , A/bh, 
and Mj rr measuring the geometrical quantities of appar- 
ent horizon. The disagreement between the event horizon 
and the apparent horizon may be large if the spacetime 
is not stationary, e.g., during the mass accretion phase in 
which the black hole mass dynamically increases. This 
fact makes the reliability of these methods worse. It 
shoul d be noted that the dynamical horizon formalism 
(e.g. JSchnetter et al.ll2006|) could be used to obtain more 
reliable estimation for mass and angular momentum of a 
dynamical black hole. 

Instead of using Eqs. (|4"5"]) and (|4*5]). we estimate an- 
gular momentum of a black hole using the conservation 
law, 



JBH = Jcon — Ju 



J, 



r>r A H 



(47) 



where Jtot is the total angular momentum of the sys- 
tem, J r >r AH is the amount of angular momentum located 
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outside the apparent horizon, and A J v is the amount of 
angular momentum carried away by neutrinos. We here 
ignore a small contribution of AJ„. Then, we adopt the 
quantity 

_ cJ con , AO \ 

= GM^ (48) 

as an approximate indicator of the non-dimensional spin 
parameter of a black hole. 

An accretion disk will be formed in the collapse of the 
rotating models. Because it is difficult to strictly define 
disk mass, we approximately estimate it by 

M disU = J p,d 3 x, (49) 

where p cu t is a cutoff density which characterizes density 
near the surface of the accretion disk, tah is radius of ap- 
parent horizon, and r cu t is a cutoff radius which charac- 
terize the size of the accretion disk. Although A/disk is no 
more than an approximate indicator, the disk mass may 
be estimated by Mdisk in a reasonable accuracy: When 
p cu t is larger than the surface density, slight change of 
p cut will result in large change of Md; s k. By contrast, 
in the case that p cu t is smaller than the surface density, 
Mdisk will not change much even if p cut is decreased to 
some extent, because density outside the disk is low. We 
choose p cut so that A/disk is not largely affected by a small 
change in p cu t and typically set p cu t = 10 10 g/cm 3 . 

In this paper, we basically consider two rates, mass 
accretion rate into a black hole (Mbh) and mass infalling 
rate onto an accretion disk (Mdisk), which are associated 
with time evolution of Mbh and Mdisk, respectively. The 
total mass infalling rate onto the system of a black hole 
surrounded by an accretion rate is then approximately 
given by M = M B h + Mdisk- 

2.6. Grid Setting 

In numerical simulations, we adopt a nonuniform grid, 
in which the grid spacing is increased according to the 
rule 

dxj+i = (1 + S)dxj, dzi + \ = (I + S)dzi, (50) 

where dxj = Xj+\ — Xj, dz\ = zi + \ — Zi, and 8 
is a constant. In add i tion, a regridding techniq ue 
(jShibata fc Shanirol 12001 ISekiguchi fc Shibatal l2005h is 
adopted to assign a sufficiently large number of grid 
points inside the collapsing core, saving the CPU time 
efficiently. The regridding is carried out whenever the 
characteristic radius of the collapsing star decreases by 
a factor of 2-3. At each regridding, the minimum grid 
spacing is decreased by a factor of ~ 2 and the geomet- 
rical factor S is changed slightly. 

All the quantities on the new grid are calculated us- 
ing the fifth-order Lagrange interpolation. However, for 
the fluid quantities such as p and h, the fifth-order in- 
terpolation could fail because the interpolation may give 
negative values of p and h — 1. In such cases, we adopt 
the linear interpolation to calculate the quantities on 
the new grid, b a sed o n the prescription proposed by 
lYamamoto et al.l (|2008). In each regridding, we solve 
the Hamiltonian constraint equation numerically. 

To avoid discarding a large amount of the matter in the 
outer region (i.e., for approximately keeping the location 
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Fig. 5. — Time evolution of the central values of density and tem- 
perature for the spherical model. The collapsing core experiences 
weak bounce at t fa 1168 ms. We note that apparent horizon is 
formed at t fa 1193 ms. 

of outer boundary) , we also increase the grid number at 
each regridding. For the regridding, we define a relativis- 
tic gravitational potential $ c = 1 — a c ($ c > 0) where 
a c is the central value of the lapse function. Because $ c 
is approximately proportional to M / R where M and R 
are characteristic mass and radius of the core, can 
be used as a measure of the characteristic length scale of 
the stellar core for the regridding. 

To check the convergence of results, a simulation in a 
finer grid resolution is also performed. Table [T] summa- 
rizes the regridding parameters (N and L are mesh num- 
ber and computational domain) of each level of the re- 
gridding procedure for normal (upper) and higher (lower) 
resolutions. 

3. RESULTS 

3.1. Spherical model 

In this section, we describe the features of collapse dy- 
namics for the spherical model as a baseline for the rota- 
tional models described later. As in the core collapse 
of an ordinary supernova for which the central value 
of entropy per baryon is s/ks ~ 1, gravitational col- 
lapse is triggered by the electron capture and the photo- 
dissociation of heavy nuclei. Then the collapse in the 
early phase proceeds in a homologous manner. Because 
of the higher value of the entropy per baryon (s/ks = 8), 
most of heavy nuclei are resolved into heliums by the 
photo-dissociation (cf. Figure E]). As the collapse pro- 
ceeds and as a result, temperature increases, heliums are 
resolved into free nucleons (p, n) . As we shall see below, 
due to the higher entropy and the resulting difference in 
the baryon composition, the collapse dynamics in a late 
phase is different from that of an ordinary supernova 
core. 

3.1.1. Gas pressure dominated bounce 

It is known that an ordinary supernova core experi- 
ences a bounce when the central density exceeds the nu- 
clear density (p nuc ~ 2 x 10 14 g/cm 3 ) above which the 
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TABLE 1 

Summary of the regridding procedure 





$ c < 0.0125 


< $ c < 0.025 


< * c < 0.05 


<<S> C < 0.1 


$ c < 0.2 


$ c > 0.2 


Axo (km) 


10.1 


4.8 


2.2 


0.98 


0.45 


0.22 
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0.008 


0.0075 


0.007 


0.0065 


0.006 


0.0065 


N 


316 


412 


524 
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L (km) 


14600 
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Fig. 6. — Evolution path of the central values of density and 
temperature for the spherical model in p-T plane (thick black solid 
curve). The thick and thin red solid curves show the boundary at 
which the condition P e = Pgas or P r = Pgas is satisfied. The thin 
black solid curves show evolution paths with constant entropy per 
baryon for s/kg = 5, 8 and 16. The two blue dashed curves denote 
the values of (p,T) with which 56 Fe or 4 He will be half by mass 
due to the photo-dissociation. An evolution path of the central 
values of density and temperature for an ordinary supernova core 
(Sckiguchi 2010b, see the text for details) is shown together (solid 
green curve). 

pressure increases drastically due to the repulsive nu- 
clear force. In the present case, the collapse is not de- 
celerated by the nuclear force but by the thermal gas 
pressure P gas at a density far below p nuc . Such a feature 
of dy namics was already reported in the recent simula- 
tions fFrver et al.ll200lHNakazato et~aIll2007t[Suwa et al.l 



12007b ). We reconfirm this previous discovery and clar 
ify the origin of this phenomena in more detail in the 
following. 

The evolution of the central values of density and tem- 
perature for the spherical model is shown in Figure [5j 
At t w 1168 ms the core experiences a weak bounce (see 
also Figure[7]). The central density at the bounce is below 
the nuclear density (« 2 x 10 12 g/cm 3 ) and the central 
value of the temperature is « 13 MeV. At these values 
of central density and temperature, the pressure in the 
inner core is dominated by the thermal pressure of gas 
composed primarily of free nucleons and heliums. 

This situation is different from that for t < 1160 ms, 
for which the pressure in most region of the inner part is 
dominated by the degenerate pressure of relativistic elec- 
trons. Because the adiabatic index of non-relativistic gas 
is r = 5/3, which is much larger than that for relativistic 
degenerate electrons, T w 4/3, the collapse is decelerated 
due to a sudden increase of the pressure. The radial pro- 
files of temperature, density, entropy, and entropy per 



baryon at the bounce along the equator are shown in 
Figure [Jj This figure shows that the profiles do not vary 
significantly after the bounce, for 1168 < t < 1183 ms. 

The critical value of entropy per baryon for the gas- 
pressure-dominated bounce to occur may be approxi- 
mately estimated as follows. We plot paths along which 
entropy per baryon is constant in Figure |5] (see the thin 
black curves). For s/ks < 5, paths of P e = P gas and 
constant entropy do not intersect. For s/ks > 16, on 
the other hand, the gas-pressure-dominated bounce can- 
not occur because the pressure is always dominated by 
the radiation pressure of photons (see the thin red curve 
in Figure [B]). Therefore, most of the results obtained in 
this paper would be applied qualitatively to models with 
5<s/fcs<16. 

3.1.2. Shock stall and Black hole formation 

As in the case of ordinary core collapse, a shock wave 
is formed at the gas-pressure-dominated bounce, and 
then, it propagates outward (see Figure [Jj. Because this 
bounce is weak, the shock wave is stalled soon after the 
bounce, at f w 1179 ms (cf. FigureEJ. Near the stalled 
shock, a region of negative gradient of electron fraction 
{dY e /dr < 0) is formed (see the blue curve in Figure [JJ 
because neutrinos carry away the lepton number from 
the shock-heated region. It is known that such a con- 
figuration is unstable to convection. However, because 
the thermally supported hot inner core quickly (~ 10 
ms) collapses to a black hole, convection does not play 
an important role by contrast with the case of ordinary 
supernovae. 

Figure |8] plots the time evolution of black hole mass for 
the spherical model. Note that the three masses of the 
black hole (see Section 12. 5[) approximately agree with 
each other (see Figure [8]) . Apparent horizon is formed 
at t 1193 ms. After the apparent horizon formation, 
we continue th e simulation using a hydrodynamic exci- 
sion technique (lHawke et al.ll2005l ). similar to adopted in 
iSekiguchi fc Shibatal (|2007f )7 

Black hole mass at the moment of its formation is 
« 5.8M , which is much larger than the maximum mass 
of cold spherical neutron stars (M co idNS,max ~ 2.2M Q 
for Shen's EOS). This is because the maximum mass of 
a hot neutron star can be much larger than the canon- 
ical value M co idNS,max due to the higher entropy. It is 
found that approximate average value of the entropy is 
s/fce ~ 7 just before the blac k hole formation (see Fig- 
ure [7]) . iNakazato et al.l (|2007f) calculated the maximum 
mass of a hot neutron star using Shen's EOS. According 
to their result, the maximum mass is ~ 5.6M Q for an 
isentropic core of s/ks ~ 7 with Y e — 0.1, which agrees 
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Fig. 7. — Radial profiles of temperature, density, entropy per baryon, and electron fraction along the radial coordinate in the equator at 
t PS 1141, 1168 (bounce), 1179 (shock stall), and 1192 (just before the apparent horizon formation) ms. Formation of a shock for t > 1168 
ms is due to the weak bounce. 
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Fig. 8. — Time evolution of black hole mass for the spherical 
model. 

approximately with our present result. After the forma- 
tion of the black hole, its mass increases gradually as the 
accretion of the material from the outer region proceeds. 
In the first ^100 ms, the mass accretion rate into the 
black hole is M B h ~ 3OM s _1 . 

3.1.3. Neutrino luminosities 

Figure [9] plots the time evolution of neutrino luminosi- 
ties for the spherical model. Before the weak bounce, av- 
erage energy of fj, and r neutrinos is largest and electron 
neutrinos are dominantly emitted and emissivity of elec- 
tron anti-neutrinos is much smaller because electrons are 
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Fig. 9. — Time evolution of neutrino luminosities for the spherical 
model. Note that the black hole is formed at t 1193 ms. 



mildly degenerate with the electron degeneracy parame- 
ter of r/ e ~ 4(> 1) and the positron fraction, responsible 
for anti-neutrino emission, is small. Note that the tem- 
perature is relatively low as T ~ a few MeV. At leading 
order, ignoring the blocking terms due to weak degener- 
acy of neutrinos, energy emission rates associated with 
the electron capture and with the positron capture are, 
respectively, written as 



QHcxXpFM, 
Q? e c cxA n F 5 (-r/ e ). 



(51) 
(52) 
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Here, th e Fermi-Dirac inte grals are approximately given 
bv fe.g- lFuTler et alJU985h 

F B (-Ti e ) i*120e-i; (53) 

FM « | + ^ + ^ e 2 + ^ - 120e-(54) 

which give, for ?7 e ~ 4, F 5 (n e ) ~ 3000 and F 5 (—rj e ) ~ 2. 
For this stage, it is found X p /X„ ~ 0.1 where X„ and 
X p are the neutron and proton fractions. Therefore, the 
relation of > Q? e c holds. 

After the weak bounce, the degeneracy parameter be- 
comes low as r\ e ~ 1.5 because high temperature of 
T > 10 MeV is achieved. In this case, F 5 (r] e ) ~ 300 
and Fc,{~rj e ) ~ 30, and electron neutrinos and electron 
anti-neutrinos are approximately identically emitted for 
X p /A„ - 0.1 because Q£ ~ Q? e c . 

The peak luminosities of electron neutrinos (« 1.8 x 
10 54 erg/s) and anti-neutrinos (1.6 x 10 54 erg/s) are 
achieved soon after the bounce (at t « 1176 ms) be- 
cause neutrinos in the hot postshock region, where the 
density is not so large that optical depth for neutrinos is 
small, are copiously emitted. These luminosities remain 
approximately constant until black hole is formed. This 
happens due to the following competing effects; as a re- 
sult of neutrino emission, thermal energy in the neutrino 
emission region is decreased, whereas as a result of com- 
pression associated with the collapse, temperature in the 
neutrino emission region is increased. 

The peak luminosities of /i and r neutrinos, on the 
other hand, are achieved just before the black hole for- 
mation. This is because the temperature significantly 
increases (see Figure [7]) due to the adiabatic compres- 
sion to enhance the pair production channel of neutrinos. 
Note that pair processes of neutrino production depend 
strongly on the temperature as Q^ n oc T 9 . Just before 
the black hole formation, luminosities of all the species 
of neutrinos become approximately identical. This shows 
that the pair production process is dominant. 

Soon after the black hole formation at t « 1193 
ms, neutrino luminosities decrease drastically because 
the main neutrino-emission region is swallowed by the 
black hole. For the spherically symmetric case, i.e., 
in the absence of an accretion disk formation, neutrino 
luminosities damp monotonically as the density of in- 
falling material decreases. The total energies emitted 
b neutrinos over the entire time of the simulation are 
E v ,tot ~ 8.3 x 10 52 , 5.2 x 10 52 , and 4.5 x 10 52 erg for 
electron neutrinos, electron anti-neutrinos, and total of 
fi and r neutrinos, respectively. 

Before closing this subsection, we briefly compare 
our results for the spherical model with those in 
iNakazato et al.l (|2007|) . who performed spherically sym- 
metric general relativistic simulations in which the Boltz- 
mann equation is solved for neutrino transfer with rele- 
vant weak interaction processes. Note that the evolution 
after the black hole formation was not followed in their 
simulations becaus e they adopted the so -called Misner- 
Sharp coordinates ([Misner fc Shard lT9645. by which evo- 
lution of black hole cannot be followed. According to 
their results for a model with the initial entropy of 
s/Ub = 7.5, the maximum neutrino luminosities achieved 
are L Ve w L 9e w 8 x 10 53 erg/s and L Vg> » 4 x 10 53 



erg/s, which are by a factor of 2-3 smaller than those 
in our results. The primary reason for this is that their 
computation was finished before the peak luminosity is 
reached due to the choice of their time coordinate, which 
is not suitable for following black hole evolution. How- 
ever, the qualitative feature of luminosity curves for each 
species of neutr i nos in our simulation agrees with that in 
INakazato et al.l (|2007l ) for the phase before the black hole 
formation. 

3.2. Moderately rotating model 

The basic features of rotational core collapse until the 
black hole formation are qualitatively the same as those 
of the spherical model: Gravitational collapse is triggered 
primarily by the photo-dissociation of heavy nuclei; the 
gas-pressure-dominated bounce occurs at a subnuclear 
density; a weak shock wave is formed at the bounce and 
is stalled quickly; a black hole is formed soon after the 
bounce in rs 30-50 ms. After the black hole formation, 
on the other hand, the dynamics of inf ailing material is 
modified by the centrifugal force; an accretion disk is 
formed around the black hole as the material with suf- 
ficient specific angular momentum falls into the central 
region. We first describe the feature of the collapse for 
the moderately rotating model in Sections 13.2.11 13.2.21 
and !3.2.51 Then, we discuss dependence of the dynamics 
of the accretion disk formation and properties of the disk 
on the amount of rotation in Section [3.4l It is found that 
the process of the accretion disk formation and properties 
of the disk depend sensitively on the amount of rotation 
initially given. 

3.2.1. Black hole and thin accretion disk formation 

In this subsection, we describe features of dynamics of 
the first ~ 200 ms after the black hole formation. We 
note that time duration of this phase depends on the 
grid resolution but the evolution process does not de- 
pend qualitatively on it. Figure [10] plots contours of den- 
sity, electron fraction, entropy per baryon, and tempera- 
ture at selected time slices around black hole formation 
epoch. As in the collapse for the spherical model, the 
weak bounce occurs at t « 1339 ms, and then, convec- 
tively unstable regions with negative gradients of electron 
fraction appear when the shock wave is stalled. However, 
because the core immediately collapses to a black hole, 
the convection is only weakly activated and plays a mi- 
nor role (see the left and middle panels in Figure fTO)) . 
Accompanied with the black hole formation, a geomet- 
rically thin, 'sub-Keplerian' disk is formed around the 
black hole (see below for details). Note that the disk is 
geometrically thin not due to the neutrino cooling (be- 
cause the disk is optically thick), but mainly due to the 
ram pressure of infalling material (see Eq. [57] and the 
discussion below) . 

Figure [TT] plots the time evolution of mass and spin 
parameter of the black hole as well as disk mass (Mdisk)- 
At t rj 1373 ms, a black hole of A/bh ~ 6.5M Q with spin 
parameter of (?bh ~ 0.6 is formed. The initial mass of the 
black hole is larger than that in the spherical collapse be- 
cause the threshold mass for the black hole formation is 
larger due to effect of the rotation (the centrifugal force). 
Note that M ce seems to be a good indicator of mass of a 
black hole even in the presence of a massive accretion disk 
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Fig. 10. — Contours of rest mass density (panels in the first row), electron fraction (panels in the second row), entropy per baryon (panels 
in the third row), and temperature (panels in the fourth row) at t as 1367 (left panels), 1374 (middle panels), and 1444 (right panels) for 
the moderately rotating model. The black regions in the contours of rest mass density and entropy per baryon, and the white regions in 
the contours of electron fraction at t = 1374 and 1444 ms are inside the apparent horizon. 



as suggested in Shibata (2007), because the time evolu- 
tion of M con and M ce approximately agrees with each 
other. The upper panel in Figure [T2l plots the time evo- 
lution of mass accretion rate into the black hole (Mbh)- 
The mass accretion rate soon (10 ms) after the black hole 

formation is high as Mbh * 4OM s _1 . The mass ac- 
cretion rate decreases gradually with time, but even at 
t - 1800 ms, it is still as high as M B h ~ 5-10M o s" 1 
(see the upper panel in Figure II 2j) . 

FigurerT3]plots the time evolution of neutrino luminosi- 
ties for the moderately rotating model. As in the spher- 
ical model, electron neutrinos are dominantly emitted 



before the weak bounce, and electron neutrinos and elec- 
tron anti-neutrinos are approximately identically emit- 
ted after the bounce. The luminosity curves of elec- 
tron neutrinos (« 1.6 x 10 54 erg/s) and anti-neutrinos 
(1.4 x 10 54 erg/s) achieve the first peak soon after the 
weak bounce (at t « 1330 ms). By contrast with the 
spherical model, the second peak appears in the neutrino 
luminosity curves at t « 1360 ms. Because oblate (or 
torus-like) neutrino 'sphere' is formed after the bounce 
due to the rotation, the optical depth of neutrinos is 
smaller in the z-direction (see the middle panels in Figure 
riTJ|) . As a result, neutrinos are more efficiently emitted 
in the z-direction and this effect constitutes the second 
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Fig. 11. — Time evolution of mass (the top panel) and the non- 
dimensional spin parameter (the lower panel) of the black hole and 
disk mass (the bottom panel) for the moderately rotating model. 
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Fig. 12. — Mass accretion rate into the black hole dM^^/dt = 
Mbh {the upper panel) and an efficiency of neutrino emission 
L v I Mbhc 2 (the lower panel) as functions of time after the black 
hole (BH) formation for the moderately rotating model. 
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Fig. 13. — Time evolution of neutrino luminosities for the moder- 
ately rotating model for the lower (solid curves) and finer (dashed 
curves) resolutions. A black hole is formed at t fa 1373 ms. 



peak. In this phase, more electron anti-neutrinos are 
emitted than electron neutrinos < Qp^) because 

the electrons inside the torus is only weakly degenerate 
7/ e ~ 1 due to high temperature and the fraction of neu- 
trons is larger than that of protons as X p /X n ~ 0.2, 
enhancing the reaction of n + e + — > p + v e . 

Soon after the black hole is formed, most of the ma- 
terial inside the oblate structure is quickly swallowed by 
the black hole because they do not have enough angu- 
lar momentum to retain in the orbit around the formed 
black hole. However, a small amount of the material with 
sufficient angular momentum forms a geometrically thin 
accretion disk around the black hole (see the right panels 
in Figure ITTJ1) . Mass of the geometrically thin disk just 
after the black hole formation is Mdisk ~ O.2A/ and sub- 
sequently decreases to be « 0.1M Q (see the bottom panel 
in Figure 1 1 If) because material with high density located 
near the rotational axis, which does not have sufficient 
angular momentum and does not constitute the disk, is 
swallowed by the black hole. Then the thin-disk mass 
relaxes to a quasi-stationary value of ~ O.IMq and the 
net mass infall rate onto the thin disk vanishes approx- 
imately (Mdi s k ~ 0). (For sudden increase of Mdi s k at 
t w 1580 ms, see Section EHH) 

The rest mass density and temperature of the thin disk 
are initially ~ 10 11 g/cm 3 and ~ 8 MeV (see Figure [TO)), 
and accordingly, the thin disk is optically thick to neu- 
trinos with the maximum optical depth of t„ ~ 4 (which 
increases as the material with high angular momentum 
falls onto the thin disk). At the same time, shocks are 
formed in the inner part of the thin disk, converting ki- 
neti c energy of infalling materi als into thermal energy 
(see ILee fc Ramirez-Ruizl (|2006f ) for discussion of a sim- 
ilar phenomenon). 

The shock is successively formed due to infall of the 
material with angular momentum not large enough to 
retain the orbit around the black hole. After hitting 
the surface in the inner region of the disk, such material 
falls into the black hole quickly because of its insufficient 
specific angular momentum, and contributes to a rapid 
growth of the black hole. A part of the thermal energy 
generated at the shock is advected together into the black 
hole (see discussion below). 

The thermal energy is also carried away by neutri- 
nos because the cooling timescale of neutrino emission, 
^cooi, is short due to the low density and small pressure 
scale height of the disk, H (although the optical depth 
is greater than unity), as 
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This is much shorter than the advection time scale ap- 
proximately given by 
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(56) 



where i?disk(~ Hsco 3* H) and v^dv are characteristic 
radius of the disk and characteristic advection velocity. 

The pressure scale height may be approximately 
determined by the following force balance relation 
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Fig. 14.— Contours of rest mass density at t m 1578 (top left), 1584 (top middle), 1591 (top right), 1644 (bottom left), 1706 (bottom 
middle), and 1800 ms (bottom right) for the moderately rotating model. 




Fig. 15. — Contours of entropy per baryon for the moderately rotating model. The selected time slices are the same as those in Figure IT4l 
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Fig. 16. — Profiles of Q/Qjf^, density, entropy per baryon, and total lepton fraction along the radial direction in the equator at t ra 1574, 
1576, 1578, and 1580 ms. 
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Fig. 17. — Contours of the total neutrino emissivity for the moderately rotating model. The selected time slices are the same as those in 
Figure [Til 
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Fig. 18. — Contours of electron fraction with velocity fields at t fa 1589 (top left panel), 1590 (top right panel), 1596 (bottom left panel), 
and 1644 ms (bottom right panel). 

where p{ and v { ~ (2GM B H/-Rdisk) 1/2 ~ 0.4-0. 5c are 
the density and velocity of the infalling material, respec- 
tively. Since -Pdisk ~ -Pram I *C Pdisk > the pressure scale 
height is very small as H/R^isk <C 1 in the early stage of 
the thin disk. 

The lower panel in Figure [T^] plots an efficiency of neu- 
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where pdisk and Pdisk are characteristic density and pres- 
sure of the disk, and P ram is the ram pressure of the 
infalling material, respectively. Equation (|57j) gives 
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1/2 
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disk 



Pdisk 



-Rdisk \ 10 30 dyn/cm 



10 10 g/crn 



(58) 



Because the density and temperature retain to be low 
due to rapid advection and copious neutrino emission, 
Pdisk ~ 10 30 dyn/cm 2 is as small as the ram pressure, 
approximately written as 



trino emission defined by L^ tot /(M-Qnc ), where L„ tot 
is the total neutrino luminosity. The efficiency is low as 
~ 0.01 in the thin accretion disk phase (until w 200 ms 
after the black hole formation). On the other hand, or- 
der of magnitude of the thermal energy generated at the 
shock in the inner region of the thin disk is estimated to 
give 
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Fig. 19 . — Contours of the Solberg-Hoiland frequency for the moderately rotating model. The selected time slices are the same as those 
in Figure Il4l 

Because the total mass of the material surrounding 
the black hole is much larger than that in the spheri- 
cal model, the neutrino luminosity remains high > 10 53 
erg/s even after the black hole formation (see Figure IT5]) . 
For ~ 200 ms after the thin disk formation, the luminos- 
ity slightly increases but is kept to be ^ 2 x 10 53 erg/s. 
Because the duration of the neutrino emission from the 
thin disk is much longer than that before the black hole 
formation, neutrinos are likely to be primarily emitted 
from the accretion disk (torus), not during the black hole 
formation, in the moderately rotating model. 
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Fig. 20. — Profiles of M v is along the radial direction in the equa- 
tor in the geometrically-thin-disk phase (at t as 1556 ms) and the 
convective phase (at t m 1644 and 1574 ms). 

where r ~ 0.1GMbh/c 2 is the distance from the black 
hole. Here, it is assumed that most of the material falling 
onto the system experiences the shock heating (i.e., the 
total mass accretion rate M is used) and an approxima- 
tion of M — Mbh + A/di s k ~ A/bh is used. 

Thus, the neutrino luminosity is by one order of magni- 
tude smaller than that the energy generated at the shock. 
This indicates that the amount of material which experi- 
ences the shock heating is much smaller than that swal- 
lowed into the black hole because of small geometrical 
cross section with the disk. As the material with high 
specific angular momentum falls onto the disk and the 
size of the disk increases, neutrino luminosities and the 
shock heating efficiency increase (see Figure [T3] and the 
lower panel in Figure [T2|) . (For decrease of neutrino lu- 
minosities at t w 1470 ms, see Section [3X2]) 



3.2.2. Disk expansion and torus formation 

Figures IT41 and [T5l plot contours of density and entropy 
per baryon at selected time slices ~ 200-400 ms after the 
black hole formation. It is found that the geometrically 
thin accretion disk formed in the early stage expands to 
form a geometrically thick accretion torus. Note that 
the disk is also 'sub-Keplerain' in this stage (see the top- 
left panel in Figure IT51) . The feature of dynamics can be 
explained as follows. 

As the material with higher specific angular momen- 
tum in the outer region falls onto the disk, the density 
and mass of the disk increases (see the bottom panel in 
Figure ITT]) . This situation is different from that in the 
early evolution of the geometrically thin disk, in which 
the material with small specific angular momentum dom- 
inantly falls. As a result, neutrino optical depth increases 
and neutrino cooling timescale becomes longer (cf. Eq. 
(|55|) V This helps further storing thermal energy inside 
the disk and the pressure scale height increases (see the 
top-left panel in Figure [T5|) . 

As the thermal energy is stored, the disk height H in- 
creases according to Eq. (|57|) . The density and the tem- 
perature (Tdisk) inside the disk eventually increase to be 
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> 10 11 g/cm 3 and > 10 MeV (and hence, P d isk > 10 30 
dyn/cm 2 ). At the same time, the ram pressure decreases 
to be < O.lPdisk (<C Pdisk) because the density of the in- 
falling material decreases to < 10 9 g/cm 3 . Consequently, 
H increases to be ~ -Rdisk (see the top-middle panels in 
Figures ll"4l and IT5|) . For H > -Rdisk > the approximate 
force balance relation ([57)) changes to 
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Because the binding due to the gravitational force by the 
black hole decreases as H increases, the disk expands 
forming a shock wave once the cond ition H > i?disk is 
satisfied (Sekiguchi & Shibata 2007). Figure [TC] shows 
that the shock is formed at t « 1576 ms. 

The neutrino opacities decrease as the disk expands 
(density and temperature decrease), and accordingly, the 
cooling timescale becomes shorter. Then, the shock 
wave is stalled and the disk relaxes to a new geomet- 
rically thick state. The shock becomes a standing ac- 
cretion shock and expands gradually because the mate- 
rial with higher specific angular momentum continuously 
falls onto the shock and also because the ram pressure 
of the infalling material continues to decrease (see the 
bottom panels in Figures [14] and [T5j • 

Note that when the pressure scale height, and thus, 
the optical depth become sufficiently large, the neutrino- 
cooling timescale becomes longer than the advection 
timescale into black hole, and consequently, neutrinos 
are trapped in the accretion flow. This can be seen 
in the time evolution of neutrino luminosities plotted 
in Figure 1131 At t ~ 1490 ms, neutrino luminosi- 
ties start decreasing slightly. The trapping of neutri- 
nos are als o found in a steady high-density accretion 
disk m odel pi Matteo et al.ll2002t iChen fc Beloborodovl 
|2007f ). Note also that similar decrease of neutrino lu- 
minosities has been found in the simulations of ordinary 
core collapse soon after the o nset of neutrino trapping 
(e.g., ILiebendorfer et al.|[200T[ ). 

Figure [T7] plots contours of the total neutrino emissiv- 
ity at selected time slices ~ 200-400 ms after the black 
hole formation. Neutrino luminosities are significantly 
enhanced after the thick torus formation. The reason for 
this is mainly that the amount of material which experi- 
ences the shock heating increases. The disk is optically 
thick to neutrinos at first and becomes optically thin as 
the disk expands. Then, neutrinos trapped inside the 
torus are emitted. This feature is somewhat similar to 
the so-called 'neutrino burst' associated with the early 
shock formation in the ordinary supernova explosion. 

After the expansion, the total luminosity reaches 
?» 2 x 10 54 erg/s because amount of material which 
experiences the shock heating significantly increases. 
Then the efficiency of neutrino emission is as high as 
Lv,tot/(MBnc 2 ) ~ 0.1 (see the lower panel in Figure 
I12[) . These agree approximately with the generation rate 
of thermal energy by infalling material on the standing 
shock, 
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where a characteristic value of M = Mbh + ^disk ~ 
10M o s _1 is adopted (see the bottom panel in Figure ITTI 
and the upper panel in Figure 1121) . The high efficiency 
indicates that neutrino optical depth is not very high for 
the neutrino-emission region and that advection of the 
thermal energy into the black hole is not very large in 
this phase because of the quick neutrino emission. 

3.2.3. Convective activities 

After the formation of the geometrically thick torus, 
convective motions are excited near the shocked region 
in the torus. The origin of the convection is explained as 
follows. 

The shock heating is more efficient in an inner part of 
the torus because kinetic energy of infalling material is 
larger (see the top-left panel in Figure [TBI . On the other 
hand, the neutrino cooling is less efficient in the inner 
part of the torus because of its higher density and result- 
ing larger optical depth. Then, the entropy per baryon 
becomes higher in the shocked inner region of the torus 
(see Figure [TS"]). and consequently, regions of negative en- 
tropy gradient along the radial direction near the equa- 
torial plane are developed. Also, because neutrinos are 
trapped and /3-equilibrium is achieved in the inner part 
of the torus, the total lepton fraction increases inward. 
These tendencies are enhanced as the accretion of the 
material with higher angular momentum proceeds. 

The condition for convective instabilities to occur is 
given by the so-called Solberg-Hoiland criterion (e.g., 
ITassoulTl978h . 

N s l = N B 2 V + k 2 < 0, (63) 

where Arv is the Brunt- Va isala frequency given by (e.g., 
ILattimer fc Mazureklll98lD 
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and k is the epicyclic frequenc y which may be written for 
nearly circular orbits as fe.g.. iBinnev fc Tremainel[l987t ) 
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Figure [Jj)] plots the profiles of angular velocity, total 
lepton fraction, and entropy per baryon along the radial 
direction in the equator after the convection sets in. It 
is clearly shown that negative entropy gradient is formed 
in several regions inside the torus, and drives convection 
(see Figures [14] and [15]) . Rotation does not play an im- 
portant role in suppressing the convective activities be- 
cause the angular velocity f2 is smaller than the Kepler 
angular velocity given by 
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(see the top-left panel in Figure ITS"]) . and thus, Coriolis 
force is not large enough. 

The convective flows cannot move freely because the 
material infalling from the outside of the torus prevents 
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the free expansion of the convective components (see the 
top- middle panel in Figure [T4J). Figure [TBI plots contours 
of electron fraction with velocity fields. Interacting with 
the thin accretion flows, a part of the convective flows 
is swerved to form finger-like structure (see the top-right 
panel in Figure II8|) . Then, the convective components 
form a swirl. Note that regions with velocity shear ap- 
pear at the interface between the convective fingers and 
the accretion flows (see the right panel in Figure [TBI . and 
hence, the Kelvin-Helmholtz instability could be devel- 
oped at the interface, generating turbulent motions (see 
the bottom- left panel in Figure [T8]) . 

In addition, oscillations of the standing shock wave are 
induced. Such shock oscillations are proposed in a dif- 
ferent contex t to explain quasi-pe riodic oscillations of X- 
ray binaries (jMolteni et al.lfl996| ) and found in a recent 
Newtonian simulatio n of sub-Kepler ian accretion flows 
around a black hole (|Giri et al.ll2010f ). 

Associated with the convective motions, many shock 
waves are formed and accretion flows show very com- 
plicated features. Because of interplay of the neutrino- 
trapping, the Kelvin-Helmholtz instability, and the con- 
vective shock, the accretion flow remains convectively un- 
stable. Figure shows the Solberg-Hoiland frequency, 
JVsh defined in Eq. (|63l) . The effective gravity appeared 
in Eq. is approximately evaluated using the New- 

tonian gravity as g c s = GM^n/r 2 . As this figure shows, 
several regions inside the standing shock remain convec- 
tively unstable. 

As a natural consequence of the convective activities 
of the accretion flow, neutrino luminosities vary violently 
in time (see Figure [T"3|) . If GRBs are driven by the pair 
annihilation of neutrinos and anti-neutrinos, such time- 
variability may explain the observed time-variability of 
GRB light curves. Furthermore, electrons in the convec- 
tive regions are only weakly degenerate due to the high 
entropy and temperature. Consequently, the emissivi- 
ties of electron neutrinos and electron anti-neutrinos are 
approximately identical (Q° c ~ Qp C ). This is favorable 
for the pair annihilation of neutrinos to electron-positron 
pairs because its rate is proportional to L U L P (see Sec- 
tion [473]). We finally note that the total energies emitted 
in neutrinos over the entire time of the simulations are 
E„,tot « 3.8 x 10 53 , 3.9 x I0 53 , and 9.4 x 10 52 erg for 
electron neutrinos, electron anti-neutrinos, and total of 
(i and r neutrinos. 

3.2.4. Effect of viscosity and formation of viscous accretion 

disk 

Finally we remark possible effects of viscosity in the 
evolution of the accretion disk, which are not taken into 
account in our simulation. Assuming that the disk (or 
torus) can b e described by the standard disk model with 
a- viscosity ()Shakura fc Sunvaevl 119731 ) . mass accretion 
rate of disk material into the black hole due to the vis- 
cous transport of angular momentum (M v ; s ) is written 
as 



M v 



(67) 



where a V i S is the viscous parameter and the pressure scale 
height is approximately estimated by 



Figure [30] plots characteristic values of M v - m along the ra- 
dial direction in the equatorial plane in the geometrically- 
thin-disk phase (at t s=s 1556 ms), early (at t ~ 1644 ms) 
and late (at t « 1772 ms) stages of the convective phase. 
During the evolution of the accretion disk, viscosity is not 
likely to play an active role as described in the following. 

In the geometrically-thin-disk phase, the predicted vis- 
cous mass accretion rate is small as M v is ^ 0.1 M Q s _1 
for a relatively large viscous parameter of a V j S = 0.1. 
Then characteristic timescale for viscous mass accre- 
tion is ~ 1 s because the disk mass is A/disk ~ O.1M 
(see Figure lll[) . which is much longer than the dura- 
tion of the geometrically-thin-disk phase ~ 200 ms.. 
Thus, the viscosity will not play an important role in 
the geometrically-thin-disk phase. 

In the convective phase, the viscous mass accretion 
rate becomes large as Af v j s ~ Af Q s" 1 for a v i s = 0.1. 
On the other hand, the mass infall rate onto the torus is 
A/disk ~ 3-4M Q s _1 (see Figure fTT|) . which is larger than 
the viscous mass accretion rate. Thus, effect of viscos- 
ity is not likely to play a central role and the disk will 
accumulate mass even in the presence of the viscosity. 

The disk will spread outward with accumulat- 
ing mass until the viscous mass accretion rate ex- 
ceeds the infall mass accretion rate onto the disk 
(A/disk ~ 47ri? 2 !isk( Ofi;f). When Mdi s k becomes smaller 
and the torus becomes more massive due to accre- 
tion of material from outer regions, the viscosity will 
play an important role on evolution and dynamics 
of the torus. Over the past decade, many groups 
have studied prope rties of the viscous accretion disk 
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successfully explained the energetics of LGRBs. 

It should be note that in the viscous accretion phase, 
the material with low angular momentum will also fall 
in the vicinity of the black hole and shock dissipation of 
the infall kinetic energy will also occur. Material with 
high angular momentum can dissipate their infall kinetic 
energy on the standing shock before they reach the cen- 
trifugal barrier. The amount of such materials depends 
on the initial density and rotational profile yet poorly 
known. There might be substantial amount of mass ac- 
cretion and energy generation due to such processes. 

3.3. Dependence on grid resolution and numerical 
accuracy 

Because the present simulation is long-term one, we 
here describe dependence of results on the grid resolu- 
tion and numerical accuracy. In Figure I13[ we compare 
the time evolution of neutrino luminosities derived both 
in the high (dashed curves) and low (solid curves) res- 
olution runs. The neutrino luminosities in the two grid 
resolutions agree very well until the black hole forma- 
tion, indicating that converged results are obtained for 
such phase. In the geometrically thin disk phase, on the 
other hand, the luminosities in the finer resolution are 
systematically higher than those in the lower resolution. 
This is because the vertical structure of the geometrically 
thin disk and shock-heated region are more accurately re- 
solved in the finer resolution, and hence, the maximum 
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Fig. 21. — Time evolution of the total baryon mass (the top 
panel), the total ADM mass (the middle panel), and the total 
angular momentum (the bottom panel) for the moderately rotating 
model. The red and blue curves correspond to the results in the 
lower resolution and in the higher resolution, respectively. 
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Fig. 22. — Time evolution (normalized) in each regrid level of the 
total baryon mass (the top panel), the total ADM mass (the middle 
panel), and the total angular momentum (the bottom panel) for 
the moderately rotating model. The solid curves correspond to the 
results in the lower resolution and the dashed curves to those in 
the higher resolution. 

temperature is higher in the finer resolution. Also, the 
geometrically thin disk more quickly expands to be the 
geometrically thick disk. This is because the thermal 
energy is more efficiently stored in the disk because neu- 
trino opacities are larger due to the higher density and 
temperature. These results indicates the importance of 
resolving the vertical structure of the geometrically thin 
disk for the quantitative study. If the grid resolution is 
not sufficient, a geometrically thin disk may remain thin 
instead of expanding to be thick torus. 

Note that the effects of grid resolution works in a pos- 
itive manner in our results, that is, the transition of a 
thin disk to a thick disk is more likely to occur. We 
therefore safely conclude that qualitative feature of our 
results does not depend on the grid resolution. 




200 400 600 800 1000 1200 1400 1600 1800 
Time [ms] 



Fig. 23. — Time evolution of the Hamiltonian constraint error 
for the moderately rotating model. The red and blue curves cor- 
respond to the results in the lower resolution and in the higher 
resolution, respectively. 
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Fig. 24. — Time evolution (normalized) in each regrid level of LI 
norm of the Hamiltonian constraint for the moderately rotating 
model. The solid curves correspond to the results in the lower 
resolution and the dashed curves to those in the higher resolution. 



To check the accuracy of our results, conservations of 
the baryon mass (M»), the ADM mass (Madm) (e.g, 
I York I [19791 ) . and the total angular momentum (J), as 
well as violations of the Hamiltonian constraint are mon- 
itored during the simulation. Figure displays the time 
evolution of these quantities. The several discontinuous 
changes correspond to the regridding procedures where 
outer low density region which does not affect the evolu- 
tion of the central region is discarded. In each regridding 
level, M*, A/adm, and J conserve well. To see this more 
quantitatively, we display the time evolution of error in 
each level of the regridding until the black hole formation 
in Figure 1211 The error is given by 
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where Qregridi denotes the conserved quantities AT*, 
■Madm j an d J in the i-th regrid level. For the purpose of 
facilitating visualization, the time is normalized by the 
duration of each regridding level. 

The error of conservation of total baryon mass grows 
monotonically in time, while it is small asO(10 -3 ). A 
part of the error is caused by the outer boundary condi- 
tions for fluid quantities where a simple copy is imposed. 



22 



The error of the ADM mass shows an oscillating behav- 
ior caused by the regridding procedure, and also is small 
as < 1%. The error in total angular momentum is also 
small as a few percent, indicating good accuracy of con- 
servation. Note that after the black hole formation, we 
start to adopt the excision procedure in solving hydrody- 
namic equations, and consequently, these quantities do 
not conserve. 

Figure [23] plots the time evolution of the Hamiltonian 
constraint error defined by iShibatal (|2003a|) 
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where where ip = e^, and A denotes the Laplacian with 
respect to jij . Namely, we use as a weight factor for 
the average. This weight factor is introduced to monitor 
whether the main bodies of the system (inner cores and 
dense matter regions), in which we are interested, are 
accurately computed or not. 

The several distinct spikes correspond to the regridding 
procedures where the Hamiltonian constraint equation is 
solved numerically. Until the black hole is formed, the 
constraint violation is very small as < 10 -2 and no signal 
of the increase is seen. After the black hole formation, 
degree of the violation becomes worse because of the ex- 
cision procedure. However, the violation is still small as 
~ 10 _1 , indicating the good accuracy of the simulation. 
Note that the integration in Eq. ([TO]) includes the inside 
the black hole. Figure [53] plots the time evolution (nor- 
malized) of the LI norm of the Hamiltonian constraint 
in each regrid level. Again, the violation does not show 
the signal of rapid increase. 

3.4. Dependence on rotation 

In this section, we describe dependence of the forma- 
tion process of the black hole and surrounding accretion 
disk, the convective activities inside the disk, and the 
emissivity of neutrinos, on the degree of initial rotation. 

3.4.1. Slowly rotating model 

In the slowly rotating model, a black hole with Mbh ~ 
6.3M© and qbh ~ 0.53 is formed at f w 1298 ms. The 
mass and spin parameter are only slightly smaller than 
those in the moderately rotating model. Figure [231 plots 
the time evolution of mass and spin parameter of the 
black hole as well as disk mass. The mass accretion rate 
into the black hole soon (10 ms) after the black hole 

formation is Mbh ~ 45M© s _1 (see the upper panel 
in Figure |2"T)) , which is slightly larger than that in the 
moderately rotating model. The spin parameter remains 
modest but gradually increases as in the moderately ro- 
tating model. 

As in the collapse of the moderately rotating model, 
a geometrically thin (but optically thick) accretion disk 
is formed soon after the black hole formation. In this 
case, a fraction of the material that forms the disk is 
smaller than that for the moderately rotating model due 
to lower specific angular momentum of fluid elements in 
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Fig. 25. — Time evolution of mass (the top panel) and the non- 
dimensional spin parameter (the middle panel) of the black hole 
and disk mass (the bottom panel) for the slowly rotating model. 
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Fig. 26. — Time evolution of neutrino luminosities for the slowly 
rotating model. 
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Fig. 27. — Mass accretion rate into the black hole dM^n/dt = 
Mbh (th e upper panel) and efficiency of neutrino emission 
L„/Mbhc 2 (the lower panel) as functions of time after the black 
hole (BH) formation for the slowly rotating model. 
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the slowly rotating model, and hence, the disk mass is 
smaller as Mdisk ~ 0.05M Q than that in the moderately 
rotating model and A/disk < (see the bottom panel 
in Figure [251) . However, a part of the material that falls 
onto the disk still produces shock waves in the inner part 
of the disk. Thermal energy generated at the shock is not 
efficiently stored in the disk in the early stage because 
most of the shocked material is advected into the black 
hole and neutrinos carry away thermal energy. Then, the 
disk remains geometrically thin for a long time (at least 
> 100 ms) after the formation of the black hole. 

Figure 1261 plots the time evolution of neutrino lumi- 
nosities. Before the black hole formation, the luminosity 
curves are similar to those in the moderately rotating 
model. It is found that the geometrically thin accretion 
disk emits ~ 10 53 erg/s by neutrinos. This magnitude 
is by factor of ~ 2 smaller than that for the moderately 
rotating model. The efficiency of neutrino emission is 
^,tot/(M B HC 2 ) ~ 0.002-0.003, which is by factor of - 3 
smaller than that for the moderately rotating model (see 
the lower panel in Figure[27|), indicating that less amount 
of material experiences the shock heating, and that more 
thermal energy is advected into the black hole before re- 
leased by neutrinos due to slower rotation and resulting 
shorter advection timescale. 

We do not find any enhancement of neutrino luminos- 
ity after the black hole formation in our simulation time. 
However, after the free-fall timescale of ~seconds, the 
material with higher specific angular momentum may 
eventually form a dense disk. Then, thermal energy may 
be stored inside the disk, and the disk may expand to 
be a geometrically thick torus when the ram pressure of 
the infalling material becomes sufficiently small. Fur- 
thermore, provided that the total mass accretion rate 
is sufficiently high as M > M e s _1 , neutrinos will be 
trapped in the inner region of the disk, and convective 
activities may set in as in the moderately rotation model 
(see discussion in Section l4~Tj) . If so, it is expected that 
neutrino luminosities are enhanced and show rapid time- 
variability as in the moderately rotating model. 

3.4.2. Rapidly rotating model 

In the rapidly rotating model, a black hole is first 
formed at t ps 1494 ms with mass of ps 6.8M q and the 
non-dimensional spin parameter of ~ 0.8. Figure [28l plots 
the time evolution of mass and spin parameter of the 
black hole together with disk mass. The spin parame- 
ter is much larger than that in the moderately rotating 
model as expected form Figure 21 

In the rapidly rotating model, the disk formation pro- 
cess is qualitatively different from that in the mod- 
erately rotating model. Figure [29] plots contours of 
rest mass density at selected time slices. The con- 
tour curve of r Ve — 5 is shown together as an ap- 
proximate boundary of occurrence of the neutrino trap- 
ping. Inside this curve, neutrinos are trapped because 

iadv(~ -Rdisk/Vadv ~ i?disk/0.1c) ~ t C ool(~ Ht u /c) for 
-Rdisk ~ 2iJ. 

In the moderately rotating model, a geometrically thin 
accretion disk is first formed, and then, it expands to 
be a geometrically thick torus. In the rapidly rotat- 
ing model, by contrast, a geometrically thick torus is 
formed immediately after the black hole formation be- 
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Fig. 28. — Time evolution of mass (the top panel) and the non- 
dimensional spin parameter (the lower panel) of the black hole and 
disk mass (the bottom panel) for the rapidly rotating model. 

cause the pressure gradient and the angular momentum 
of the fluid near the equator are large enough that it re- 
tains an orbit outside the ISCO. The disk at this phase 
is still 'sub-Keplerian' with fi/fijf rs 0.8 at its maximum 
and the pressure gradient plays a role in the immediate 
torus formation. Reflecting the torus formation, Mdisk 
is much larger as s=s O.4M than that in the slowly and 
moderately rotating models (see the bottom panel in Fig- 
ure [28]) ■ Shock waves formed at the weak bounce are not 
swallowed into the black hole and a torus-shaped stand- 
ing accretion shock remains around the black hole. 

Associated with the torus formation, the mass accre- 
tion rate into the black hole just after the black hole for- 
mation shows non-monotonic behavior by contrast with 
the slowly and rapidly rotating models (see the upper 
panel in Figure I30[) . The mass accretion rate quickly 

drops to be Mbh ~ 20M Q s -1 at t « 6 ms after the black 
hole formation because of the centrifugal and pressure- 
supported hangup of the torus. The subsequent oscil- 
lating behavior is due to mass accretion associated with 
the oscillation of the torus. Then the mass accretion rate 
decreases quickly with time because the centrifugal force 
of the infalling material prevents the rapid accretion into 
the black hole. Note that the pressure gradient plays a 
role also in this phase. The mass accretion is expected 
to cease when Mbh ~ 12M Q . 

Figure [3T] plots the neutrino luminosities as a func- 
tion of time. Until the onset of the weak bounce (until 
the first local peak), the luminosity curves are similar to 
those in other models. After the weak bounce occurs, the 
material near the rotation axis starts collapsing, and as 
a result, the temperature increases due to compression 
and the optical depth near the rotation axis relatively 
decreases. Then, second local peak (at t ps 1475 ms) as- 
sociated with a substantial emission from the vicinity of 
rotation axis appears. This is the same feature as found 
in the slowly and moderately rotating models. In the 
rapidly rotating model, in addition, third local peak ap- 
pears just before black hole formation at t ps 1494 ms. 
This is due to the fact that a dense torus, which subse- 
quently falls into the black hole, is formed (see the first 
panel in Figure [29]) and emits a large amount of neutrinos 
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Fig. 29.— Contours of rest mass density at t w 1495 (top left), 1497 (top middle), 1499 (top right), 1500 (bottom left), 1502 (bottom 
middle), and 1535 ms (bottom right) for the rapidly rotating model. The green curves indicate the region where t„ = 5 for electron 
neutrinos. 
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Fig. 31. — Time evolution of neutrino luminosities for the rapidly 
rotating model. Note that the black hole is formed at t as 1494 ms. 



Fig. 30. — Mass accretion rate into the black hole dMBu/dt = 
Mbh (the upper panel) and efficiency of neutrino emission 
L^/Mbh c2 (the lower panel) as functions of time after the black 
hole (BH) formation for the rapidly rotating model. 

just before swallowed by the black hole. 

After the black hole formation, the luminosities de- 
crease slightly. However, a dense torus surrounding the 
black hole is formed in a short time scale. Then, the lu- 
minosity increases again, and becomes as large as the sec- 
ond and third peaks with the total luminosity ~3x 10 54 
erg/s. Approximate generation rate of thermal energy 
at the shock on the surface of the torus due to infalling 
material is 



GM W M 



'A x 10 54 
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2OM s~ 



erg/s. (72) 



Thus, the neutrinos are emitted by converting infall ki- 
netic energy of the material to the thermal energy. 

Convective motions are also observed in the rapidly 
model as in the moderately rotating model. A large- 
scale circulation is formed associated with the formation 
of the thick, (mainly) centrifugally supported torus (see 
the bottom-middle panel in Figure l29| . However, suc- 
cessive large-scale circulations, appeared in the moder- 
ately rotating model, do not occur in the rapidly rotat- 
ing model, although small-scale convective activities are 
driven (see the bottom-right panel in Figure |2T)1) . This is 
due to the stabilizing effect of the epicyclic frequency (see 
Eq. d63j)). Figure 1321 plots the Brunt- Vaisala, frequency 
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Fig. 32. — Contours of the Brunt- Vaisala frequency (left panel) and the Solberg-Hoiland frequency defined by Eq. (163 I t (right panel) at 
t sa 1534 ms for the rapidly rotating model. 



(see Eq. (|64)) ) and the Solberg-Hoiland frequency denned 
by Eq. (|63j) . As shown in this figure, there exist regions 
with negative gradients of entropy per baryon and lep- 
ton fraction (j V R y < 0) inside the thick torus (see the left 
panel in Figure [35]). However, most of the low- frequency 
modes are suppressed by the stabilizing epicyclic mode 
and only the higher- frequency modes are present. Con- 
sequently, large-scale circulation modes are suppressed 
and only small-scale convective modes appear. 

Due to the absence of large-scale convective modes, ef- 
fects of the convection on neutrino luminosities are likely 
to be minor. Indeed, no violent time-variability is ob- 
served after the thick torus formation. The small bumps 
in luminosities at t w 1500-1510 ms are associated with 
the large-scale circulation (see the bottom-middle panel 
in Figure EH). 

The total mass of the torus is ~ 7% of the black hole 
mass and gradually increases (see Figure |55J). The self- 
gravity of the torus may play a role in a later phase; the 
torus may be unstable against non-axisymmetric pertur- 
bation and this may affect evolution of the torus because 
angular momentum transport and redistribution inside 
the torus are enhanced. To strictly clarify the evolution 
of such massive torus, a three-dimensional numerical sim- 
ulation may be needed. This is one of the issues left for 
the future. 

Finally, we remark possible effects of viscosity in the 
rapidly rotating model. Assuming that the torus can be 
described by the standard disk model, the mass accre- 
tion rate associated with a hypothetical viscous stress is 
estimated as M v i s ~ 3-5M© s" 1 for a V i S = 0.1 (cf. Eq. 
([67])). Because the mass infalling rate onto the torus is 
-Mdisk » 8M & s _1 at the late phase (see the bottom panel 
in Figure [28]), the viscosity is not expected to play a cru- 
cial role for the evolution of the torus at an early phase 
simulated in this paper. 

However, in a later phase, when the mass infalling rate 
onto the torus becomes smaller, the viscosity is expected 
to play an important role. Then, an ADAF-type (ac- 
cretion dominated accretion flow) accretion flow may be 



the outcome in the presence of a large viscosity. A high- 
velocity outflow may be accompanied bec ause the accre- 
tion rate is likely to be very high (e.g., iNaravan et al.1 
120011 ). The high black hole spin may also play an impor- 
tant role for driving a high-velocity outflow because the 
heating rate is enhanced near the ISCO and the mass ac- 
cretion is suppressed due to the small black hole radius. 

4. DISCUSSIONS 

4.1. Effect of the black hole spin on disk property and 
neutrino emissivity 

The black hole formed after the core collapse is in gen- 
eral not a Schwarzschild black hole but a rotating black 
hole (<zbh ^ 0.5 for our models). In addition, a high spin 
state with (?bh <; 0.8 is easily achieved during the evolu- 
tion of the black hole. Thus, it is necessary to take into 
account the effects associated with such a high black hole 
spin to build plausible models in the collapsar scenario. 

The spin of a black hole is known to play a 
crucial role on the evolu tion of the accretion disk 
(jChen & Beloborodovll2007| ). The inner edge of the disk 
(or torus) around a rapidly rotating black hole comes 
closer to the black hole than that around a Schwarzschild 
black hole, and consequently, the temperature and den- 
sity of the disk reach higher values. These significantly 
enhance neutrino luminosities. In addition, due to the 
higher density and temperature, the disk becomes more 
opaque to neutrinos, and neutrinos are often trapped in 
the inner regions of the disk. This leads to formation of 
regions with negative entropy gradient, and convection is 
induced. As a result of convection, neutrino luminosity 
curves may become highly variable. 

Here, it should be noted that the trapping of neu- 
trinos and occurrence of convective motions are not 
likely to be special consequences of the high mass ac- 
cretion rate (M ~ 10M© s _1 ) achieved in our mod- 
els. A^TOr^ing_toj^sjolts of a g eneral relativistic study 
by iChen fc Beloborodovl ()2007[ ) the neutrino trapping 
occurs even with a moderate mass accretion rate of 
M ~ Mq s _1 for accretion flows around a rapidly ro- 
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Fig. 33. — The rip-component of the Reynolds stress tensor t rip 
normalized so that the maximum amplitude of negative sign is 
unity. 

tating Kerr black hole. For accretion flows around a 
Schwarzschild black hole, by contrast, the neutrino trap- 
ping does not occur even with a high mass accretion rate 
of M ~ 1OM s" 1 (|Chen fc Beloborodovl 121)071) . This 
illustrates that the black hole spin plays a crucial role 
on the properties of accretion flows around a black hole. 
They also find that the neutrino trapping occurs in the 
vicinity of the black hole (r < 20GA/ B h/c 2 ), as m our 
case. This indicates the importance of resolving the re- 
gions in the vicinity of the black hole because the seed of 
convection is formed there. (We note that the enhance- 
ment of neutrino luminosities due to the convection was 
not found in previous pseudo-Newtonian studies because 
a rather wide region near the black hole was excised in 
these studies.) 

4.2. Comparison with CDAFs 

Presence of convective accretion flow, named as 
convection-do minated accret i on flo w (CDAF), was first 
predicted by iNaravan fc Yil ()1994h in their studies of 
a self-similar solution of advection-dominated accretion 
flows (ADAFs). Later, CDAFs were foun d in numeri- 
cal studies of ADAFs around a black hol e ()Stone et al.l 
Il999t Hgumenshchev fe Abramowiczl 12000ft . They found 
as a remarkable property of CDAF that the convection 
transports the angular momentum inward rather than 
outward. 

To see whether this is the case in the present simu- 
lation, we calculate the rcp-component of the Reynolds 



stress tensor, t rip = {Sv r 5v v ), where dvi 



(vi) is 



the velocity fluctuation and ( ) den otes time-averaging 
ijlgumenshchev fe Abramowicd [2000). Note that nega- 
tive (positive) sign of t rip corresponds to the inward (out- 
ward) transfer of the angular momentum. Figure [3"3"l plots 
contour of t rtp in the x-z plane. This figure clearly shows 
that there are regions with negative values of t rip near 
the outer surface of the torus. Convection in these re- 
gions transports the angular momentum inward, gener- 
ating flows with higher angular momentum in an inner 
region. Such flows will then move outward forming cir- 
culations. 

While the CDAF-like accretion flows are formed in the 
outer part of the torus, flows in the inner region are 



similar to those of neutrino-d ominated accretion flows 
(NDAFs) (jPopham et al.lll999ft . Furthermore, the torus 
is accompanied by the quasi-radial flows which consist 
of the material with low angular momentum and the 
outer geome trically thin ac c retion flows near the equa- 
torial plane. INaravan et al.l ()2001h found that transition 
between CDAF and NDAF are determined by a charac- 
teristic radius r out : Flows injected from r > r out form 



CDAFs, and those injected from 



< 



r ut 



form NDAFs. 



In terms of the specific angular momentum, transition 
between CDAF and NDAF may be determined by a char- 
acteristic specific angular momentum jout^ The mate- 
rial with j > j ou t form CDAFs and those with j < j out 
form NDAFs. As found in the present simulation, the ac- 
cretion flows in the moderately rotating collapsar model 
will be characterized by the inner NDAF-like and outer 
CDAF-like parts. 

4.3. Application to Gamma-ray bursts 

We now turn to application of our results to LGRBs. 
We consider, as two possible ways of the energy depo- 
sition process, the neutrin o pair annihilation and th e 
Blandford-Znajek process (Blandfo rd fe Znajeld 11977ft . 
Because both of which processes are not included in our 
numerical simulation, we give an order estimate of the 
energy deposition rates for the purpose of clarifying the 
potential of driving relativistic jets in our models. 

The annihilation rate of neutrinos and anti- 
neutrinos into electron-positron pairs has been 
calculated as a mechanism to p ower GRBs by sev- 
eral groups (iRuffert et alj 119971 : JPopham et all 119991: 
Asano fc Fukuvamal 120001 12001b fSalmonson fc Wilson! 



20011: iSetiawan et all 12004 120061: iBirkl et al l 120071: 



Harikae et alJ l2010al lbl: iZalamea fe Beloborodovl 120111 ) 



The energy from the neutrino pair annihilation should 
be deposited in a baryon-poor region in order to generate 
highly relativistic outflows. The funnel region near the 
rotational axis above the torus is a promising place for 
this purpose. 

Here, we an order estimate of the total energy de- 
position rate by the neutrino pair annihilation (E v „). 

The deposition ra te is proportional to M 9 / 4 M B ^ 2 
(|Beloborod ov 2008). In this estimation, the neutrino lu- 
minosity is assumed to be originated from a viscous heat- 
ing. In our present calculation, the neutrino luminosity 
is determined by mass accretion rate of the infalling ma- 
terial which experiences the shock heating and increases 
thermal energy of the disk. However, the dependence of 
the pair-annihilation rate on the mass infall rate M is 
essentially the same for thick torus phase. Due to this 
strong dependence on the mass accretion rate, the energy 
deposition by the neutrino pair-annihilation will be im- 
portant only for an early phase of the LGRB formation. 

In the geometrically thin disk, the efficiency of the neu- 
trino pair annihilation for a rapidly rotating black hole 
may be written, according to a rec e nt gen eral relativistic 
study by Zala mea fc Beloborodovl (|201lft , as 



(eff)„ p = ~ 0.01 

l^v. tot. 



M 



5/4 



Mrs- 1 



/ M BH 
\WM e 



-3/2 



(73) 



where L^tot is the total neutrino luminosity. In the 
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present simulation, the expected energy deposition rate 
by neutrino pair annihilation is quite high as E v9 ~ 10 53 
erg/s for Af BH ~ 10M Q , M ~ 1OM s" 1 , and L^tot ~ 
10 erg/s in an early phase of disk evolution for ~ 1 s. 

The efficiency of the neutrino pair annihilation de- 
pends strongly on the geometry of the disk. In par- 
ticular, (eff)„p is proportional to K~ 4 , where V^nn is 
chara c teristic vo l ume above the disk (iMochkovitch et al.l 
1993t iLhi et all I2010t iZalamea fe Beloborodovi 1201 ID . 
Liu et all ( 2010H calculated the vertical structure of 
geometrically-thick accretion flows in the pseudo- 
Newtonian gravity and estimated the energy deposition 
rate. They found that the efficiency could be enhanced 
by an order of magnitude. In this case, a very large 
energy deposition rate by neutrino pair annihilation of 
E V v ~ 10 54 erg/s may be expected. 

The outgoing Poynting power at the hori- 
zon in the Blandford-Znajek process is given by 
(Bla ndford fe Znaield[l~977l iThorne et aHll986ft 



E 



BZ 
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n B (n H -n B ) 



(74) 



where Bjj is magnitude of magnetic fields normal to the 
horizon, Rh ~ GMbh/c 2 is the radius of the horizon, 
and £1 h and 51 b are the angular velocities of the horizon 
a nd the magnetic fi eld lines. 

IMcKinnevI J2005) suggested an approximate fitting for- 
mula for the estimation of the field strength based on re- 
sults of general relativistic magnetohydrodynamical sim- 
ulations. According to his formula, the outgoing Poynt- 
ing power in the Blandford-Znajek process is given by 



E B z ~ 10 fn H 9bh 



M 



M^s" 1 



erg/s, (75) 



where fn H is a parameter which depends strongly on the 
angular velocity and the most optimistic cond ition 51b = 
51g/2 is assumed. According to the result of IMcKinnevI 
(2005), > 10% of the total outgoing power may be used 
to produce the LGRB jet. Thus, the outgoing jet power 
will be ^BZjct ~ lO^fovqBniM/iMQS- 1 )) erg/s for 
our models. 

The Blandford-Znajek power will eventually become 
much larger than the deposition rate by the neu- 
trino pair-annihilation because the power depends more 
weakly on the mass accretion rate. Even in a late phase 
with M ~ 0.1M Q s -1 , a jet power of EBZjct ~ 10 51 erg/s 
may be achieved if the black hole is sufficiently rapidly 
rotating (<7bh ^ 0.9 for which fa H > 10), accumulating 
the angular momentum of infalling material. Note also 
that magnetic fields may be amplified in the torus due 
to the magneto-rotational inst ability and/or convection 
(IBalbus fe HawlevlH99lL [1991 . 

4.4. Gravitational waves from anisotropic neutrino 
emission 

A cosmological population of core-collapse supernovae 
is one of the mos t important sources o f gravitational 
wave backgrounds (jBuonanno et alj20 05) . Gravitational 
waves (GWs) associated with anisotropic neutrino emis- 
sion are particularly important because they generate a 
burst of GWs accompanying with the memory effect, 



the s o-called burst with memory (TBraginskii fe Thornd 
I1987D . GW memory due to anisotropic neutrino emis- 
sion could contaminate, at l ow frequencies around 
0.1Hz, the inflationary GW (IBuonanno et al.l 120051: 
IHiramatsu et al.l 120051: ISuwa et al.ll2007al) . which is one 
of the tar gets of future sp ace GW de tectors such as 
DEC IGO (ISeto et al.1 l200l and BBO (lUngarelli et al.l 
120051) . Here, we give an order estimate of the amplitude 
of GWs associated with anisotropic neutrino emission. 

The amplitude of GWs d ue to anisotropic neu- 
trino emission is given by (IMueller fe Jankal 119971 : 
IKotake et aL1l2007l: ISuwa fe Murasell2009fl . Taking char- 
acteristic values of total neutrino luminosity of ~ 10 54 
erg/s from our simulation results and assuming a dura- 
tion of neutrino emission of At u ~ 1 s (cf. the moderately 
rotating model), the amplitude may be estimated as 



2 x 10" 



lOGpc 
D 



10 54 erg/sy V Is 



. (76) 



where D is the distance to the source. This value is 
as large as that calculated by ISuwa et~aH (|2007aft for the 
collapse of 300Af© PopIII st ellar core colla pse. Note that 
the initial core mass in ISuwa et al.l (|2007aD is about three 
times larger than ours. The peak neutrino luminosities 
achieved in their results are by a factor of ~ 10 larger 
than those in our results, while the duration in their re- 
sults is by a factor of ~ 10 shorter than that in the 
moderately rotating model, because they failed to find 
convective activities in the accretion torus. 

If long-term neutrino emission as found in the present 
simulations is universal for the Pop III stellar collapse, 
the GW memory due to anisotropic neutrino emission 
could significantly contaminate the inflationary GW. 

5. SUMMARY 

In this paper, we performed axisymmetric simula- 
tions of very massive stellar core collapsing to a sys- 
tem composed of a rotating black hole and surround- 
ing disk in full general relativity. We took into account 
a nuclear-theory-based finite-temperature EOS (Shen's 
EOS), weak interaction processes such as electron cap- 
ture and pair-neutrino processes, and neutrino cooling, 
which is handled by a general relativistic leakage scheme 
(lSekiguchill20l7M ). 

Progenit or models of LG RBs suggested in the litera- 
tures (e.g. JFrver et alll2007l )) raise a possibility that they 
may have an entropy higher than that of ordinary super- 
nova cores. In this work, we employed a core with a high 
entropy of s/ks = 8 as the initial condition. Because the 
distribution of angular momentum in very massive stars 
is highly uncertain, we employed four models (spheri- 
cal, slowly rotating, moderately rotating, and rapidly ro- 
tating models) by superimposing a profile of rotational 
angular velocity in a parametric manner. The initial 
models adopted in this paper are not rapidly rotating 
in the sense that the rotation velocity imposed is much 
smaller than that required to retain the IS CO around a 
Schwarzschild bl ack hole and that considered in previ- 
ous studies (e.g ., MacF adven fc Wooslevi (|1999t ). see also 
lLopez-Camara et al.l ([20091 ) and references therein). 

As in the collapse of ordinary supernova cores, gravita- 
tional collapse sets in due to photo-dissociation of heavy 
nuclei and electron capture. However, the collapse dy- 
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namics and properties of neutrino emission are different 
from those of ordinary supernova cores. The characteris- 
tics of the collapse of high-entropy cores are summarized 
as follows: 

1. The gravitational contraction is decelerated by 
the thermal gas-pressure of free nucleons at a 
subnuclear density and the core experiences a 
weak bounce (the gas-pressure-dominated bounce) . 
This is a result of the high en tropy. We re- 
confirmed this previous d iscovery ()Nakazato et al.l 
120071: iSuwa et al.l l2007bD and clarified the physi- 
cal origin in detail: We clarified that the weak 
bounce is universal for the collapse of the core with 
s/k,B ~ 5-16. 

2. Because the gas-pressure-dominated bounce is too 
weak to halt the infalling material, a black hole 
is formed soon after the bounce (within ~ 30 ms). 
The mass of the black hole at the moment of its for- 
mation (~ 5.8-7M Q ) is much larger than the max- 
imum mass of a cold neutron star (w 2.2M for the 
Shen's EOS). This is also due to the high entropy 
(high thermal pressure) . Just before the black hole 
formation, the pair-neutrino production processes 
are enhanced because the temperature increases 
due to the adiabatic compression (due to neutrino 
trapping). As a result, approximately the same 
amount of electron neutrinos and anti-neutrinos are 
emitted. The mass accretion rate into the black 
hole just after the black hole formation and the to- 
tal neutrino luminosity just before the black hole 
formation are ~ 40M o s" 1 and - 4 x 10 54 erg/s 
depending weakly on the degree of rotation. Thus 
the maximum efficiency for the neutrino emission 
is L v /{Mc 2 ) ~ 6%. 

3. In the moderately rotating model, a geometrically 
thin accretion disk is first formed around the black 
hole and shocks are formed on its surface, gener- 
ated by the infalling material. As the thermal en- 
ergy is stored in the disk, it expands eventually to 
be a geometrically thick accretion torus. After the 
thick torus formation, convective a ctivities, which 
are s imilar to those in CDAFs (Nar avan et al.l 
l200l . set in because a region with negative entropy 
gradient emerges in the inner part of the torus, 
due to occurrence of the neutrino trapping. The 
neutrino luminosities are L„ c + L Pe ~ 10 54 erg/s, 
and show violent time- variability. Here we empha- 
size that the source of thermal-energy generation, 
which is eventually dissipated by neutrinos, is the 
shock heating of infalling materials. The high spin 
of a black hole is likely to play a crucial role on the 
evolution of the accretion disk, convective activi- 
ties, and the enhancement of neutrino luminosities. 

4. The evolution process of the accretion disk and 
neutrino emissivity depend strongly on the degree 
of initial rotation. In the slowly rotating model, 
the disk remains geometrically thin for a long time, 
and hence, the neutrino emissivity also remains rel- 
atively small (L ~ 10 53 erg/s) for more than 100 



ms. In the rapidly rotating model, by contrast, 
a geometrically thick torus is immediately formed 
after the black hole formation, and luminosities of 
neutrinos emitted from the torus are as high as 
10 54 erg/s even at its formation. However, the 
convection is suppressed by the stabilizing epicyclic 
mode due to the rapid rotation and no violent time- 
variability is observed in the neutrino luminosities. 

5. Irrespective of the degree of rotation, long-lived 
disk or torus surrounding the black hole is a pri- 
mary emitter of neutrinos because of its high lu- 
minosity and long lifetime J> Is. This im- 
plies that anisotropic emission of neutrinos comes 
mainly from the accretion disk (torus) surrounding 
a black hole, not from the dense matter collaps- 
ing to a black hole. For a correct estimation of 
gravitational- wave background by anisotropic neu- 
trino emission, it may be necessary to understand 
the physical condition of the accretion disk or torus 
(see below). 

Finally, we comment on major limitations of the 
present study. First, we adopt initial conditions which 
are not based on latest theoretical models of stellar evo- 
lution. We are going to perform simulations adopting 
more realistic initial models soon. Second, the present 
simulations are performed on the assumption of axial 
symmetry. The accretion disk formed in the present sim- 
ulations may beco me unstable against non - axisymmetric 
insta b ilities (e.g., Korobkin et al.l 120 lit iTavlor et al.l 
1201 U IKiuchi et al.l 20111) . Competition between non- 
axisymmetric instabilities and convective instabilities 
should be explored. Third, we do not take account of 
the neutrino heating. A simple approximated procedure 
of including effe c ts of neutrino heating is adopted by 
lO'Connor & ~Ott] (poll in which stellar core collapse to 
a black hole is studied by a spherically symmetric fully 
general relativistic simulation. We also plan to study 
effects of ne utrino heating using a recently developed 
formulation ()Shibata et al.l 1201 lh . Fourth, we do not 
consider effects of magneti c fields which will play a role 
during the collapse (e.g.. | Barkov fc Komissarovl 120081 : 



IKomissarov fe Barkovl [20091) and subsequent evolution 
of th e disk (e.g., iPenna et al.l I2010T : iBarkov fc Baushevj 
1201 If ) if progenitor cores have large magnetic fields. We 
plan to perform simulations taking account of magnetic 
fields using a general rel ativistic magnetohydrodyn amic 
code we have developed (|Shibata fc Sekiguchill2005l ). 
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